subroutine c_swap ( x, y ) ! !******************************************************************************* ! !! C_SWAP swaps two complex values. ! ! ! Modified: ! ! 26 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, complex X, Y. On output, the values of X and ! Y have been interchanged. ! implicit none ! complex x complex y complex z ! z = x x = y y = z return end subroutine cchdc ( a, lda, p, work, ipvt, job, info ) ! !******************************************************************************* ! !! CCHDC computes the Cholesky decomposition of a positive definite matrix. ! ! ! Discussion: ! ! A pivoting option allows the user to estimate the ! condition of a positive definite matrix or determine the rank ! of a positive semidefinite matrix. ! ! Parameters: ! ! on entry ! ! a complex(lda,p). ! a contains the matrix whose decomposition is to ! be computed. onlt the upper half of a need be stored. ! the lower part of the array a is not referenced. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer P, the order of the matrix. ! ! Workspace, complex WORK(P). ! ! ipvt integer(p). ! 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 ! value of ipvt(k). ! ! if ipvt(k) > 0, then x(k) is an initial ! element. ! ! if ipvt(k) == 0, then x(k) is a free element. ! ! if ipvt(k) < 0, then 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 l >= k. ipvt is not referenced if ! job == 0. ! ! job integer. ! job is an integer that initiates column pivoting. ! if job == 0, no pivoting is done. ! if job /= 0, pivoting is done. ! ! on return ! ! a a contains in its upper half the cholesky factor ! of the matrix a as it has been permuted by pivoting. ! ! ipvt ipvt(j) contains the index of the diagonal element ! of a that was moved into the j-th position, ! provided pivoting was requested. ! ! info contains the index of the last positive diagonal ! element of the cholesky factor. ! ! for positive definite matrices info = p is the normal return. ! for pivoting with 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. ! implicit none ! integer lda integer p ! complex a(lda,p) integer i integer info integer j integer job integer jp integer ipvt(p) integer jt integer k integer kb integer km1 integer l real maxdia integer maxl logical negk integer pl integer plp1 integer pu logical swapk 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 = ipvt(k) > 0 negk = ipvt(k) < 0 ipvt(k) = k if (negk) ipvt(k) = -ipvt(k) if ( swapk ) then if (k /= pl) then call cswap(pl-1,a(1,k),1,a(1,pl),1) temp = a(k,k) a(k,k) = a(pl,pl) a(pl,pl) = temp a(pl,k) = conjg(a(pl,k)) plp1 = pl + 1 do j = plp1, p if (j < k) then temp = conjg(a(pl,j)) a(pl,j) = conjg(a(j,k)) a(j,k) = temp else if (j /= k) then call c_swap ( a(k,j), a(pl,j) ) end if 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,a(1,pu),1) call c_swap ( a(k,k), a(pu,pu) ) a(k,pu) = conjg(a(k,pu)) do j = k + 1, p if ( j < pu ) then temp = conjg(a(k,j)) a(k,j) = conjg(a(j,pu)) a(j,pu) = temp else if ( j /= pu ) then call c_swap ( a(k,j), a(pu,j) ) end if end if end do call i_swap ( ipvt(k), ipvt(pu) ) 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 ( k >= pl .and. k < pu ) then do l = k+1, pu if (real(a(l,l)) > maxdia) 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.0e0 ) then info = k - 1 return end if ! ! Start the pivoting and update IPVT. ! if ( k /= maxl ) then km1 = k - 1 call cswap(k-1,a(1,k),1,a(1,maxl),1) a(maxl,maxl) = a(k,k) a(k,k) = cmplx(maxdia,0.0e0) jp = ipvt(maxl) ipvt(maxl) = ipvt(k) ipvt(k) = jp 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.0e0) a(k,k) = work(k) do j = k+1, p if ( k /= maxl ) then if (j < maxl) then temp = conjg(a(k,j)) a(k,j) = conjg(a(j,maxl)) a(j,maxl) = temp else if (j /= maxl) then call c_swap ( a(k,j), a(maxl,j) ) end if end if end if a(k,j) = a(k,j)/work(k) work(j) = conjg(a(k,j)) temp = -a(k,j) call caxpy(j-k,temp,work(k+1),1,a(k+1,j),1) end do end do return end subroutine cchdd ( r, ldr, p, x, z, ldz, nz, y, rho, c, s, info ) ! !******************************************************************************* ! !! 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 ! determineds 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. ! ! Parameters: ! ! on entry ! ! r complex(ldr,p), where ldr >= p. ! r contains the upper triangular matrix ! that is to be downdated. the part of r ! below the diagonal is not referenced. ! ! Input, integer LDR, the leading dimension of R. ! ! Input, integer P, the order of the matrix. ! ! x complex(p). ! x contains the row vector that is to ! be removed from r. x is not altered by cchdd. ! ! z complex(ldz,nz), where ldz >= p. ! z is an array of nz p-vectors which ! are to be downdated along with r. ! ! Input, integer LDZ, the leading dimension of Z. ! ! nz integer. ! nz is the number of vectors to be downdated ! nz may be zero, in which case z, y, and rho ! are not referenced. ! ! y complex(nz). ! y contains the scalars for the downdating ! of the vectors z. y is not altered by cchdd. ! ! rho real(nz). ! rho contains the norms of the residual ! vectors that are to be downdated. ! ! on return ! ! r ! z contain the downdated quantities. ! rho ! ! c real(p). ! c contains the cosines of the transforming ! rotations. ! ! s complex(p). ! s contains the sines of the transforming ! rotations. ! ! info integer. ! info is set as follows. ! ! info = 0 if the entire downdating ! was successful. ! ! info =-1 if r could not be downdated. ! in this case, all quantities ! are left unaltered. ! ! info = 1 if some rho could not be ! downdated. the offending rhos 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) complex 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)) do j = 2, p s(j) = conjg(x(j)) - cdotc(j-1,r(1,j),1,s,1) s(j) = s(j)/conjg(r(j,j)) end do norm = scnrm2(p,s,1) if ( norm >= 1.0e0 ) then info = -1 return end if alpha = sqrt(1.0e0-norm**2) ! ! Determine the transformations. ! do ii = 1, p i = p - ii + 1 scale = alpha + cabs(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 = (0.0e0,0.0e0) 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 = cabs(zeta) if ( azeta > rho(j)) then info = 1 rho(j) = -1.0e0 else rho(j) = rho(j)*sqrt(1.0e0-(azeta/rho(j))**2) end if end do return end subroutine cchex ( r, ldr, p, k, l, z, ldz, nz, c, s, job ) ! !******************************************************************************* ! !! 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 ! ! trans(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 users 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, i.e. 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. ! ! 1. right circular shift (job = 1). ! ! 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. ! ! 2. left circular shift (job = 2). ! the columns are rearranged in the following order ! ! 1,...,k-1,k+1,k+2,...,l,k,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. ! ! Parameters: ! ! on entry ! ! r complex(ldr,p), where ldr>=p. ! r contains the upper triangular factor ! that is to be updated. elements of r ! below the diagonal are not referenced. ! ! Input, integer LDR, the leading dimension of R. ! ! Input, integer P, the order of the matrix. ! ! k integer. ! k is the first column to be permuted. ! ! l integer. ! l is the last column to be permuted. ! l must be strictly greater than k. ! ! z complex(ldz,nz), where ldz>=p. ! z is an array of nz p-vectors into which the ! transformation u is multiplied. z is ! not referenced if nz = 0. ! ! Input, integer LDZ, the leading dimension of Z. ! ! nz integer. ! nz is the number of columns of the matrix z. ! ! job integer. ! job determines the type of permutation. ! job = 1 right circular shift. ! job = 2 left circular shift. ! ! on return ! ! r contains the updated factor. ! ! z contains the updated matrix z. ! ! c real(p). ! c contains the cosines of the transforming rotations. ! ! s complex(p). ! s contains the sines of the transforming rotations. ! 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 km1 integer l integer lm1 integer lmk complex r(ldr,p) complex rjp1j complex s(p) complex t complex z(ldz,nz) ! ! Initialize ! km1 = k - 1 lmk = l - k lm1 = l - 1 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, lm1 j = lm1 - jj + k do i = 1, j r(i,j+1) = r(i,j) end do r(j+1,j+1) = (0.0e0,0.0e0) end do do i = 1, km1 ii = l - i + 1 r(i,k) = s(ii) end do ! ! Calculate the rotations. ! t = s(1) do i = 1, lmk 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, lmk 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, lmk 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 = lmk + i s(ii) = r(i,k) end do do j = k, lm1 do i = 1, j r(i,j) = r(i,j+1) end do jj = j - km1 s(jj) = r(j+1,j+1) end do do i = 1, k ii = lmk + i r(i,l) = s(ii) end do do i = k+1, l r(i,l) = (0.0e0,0.0e0) end do ! ! 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, lm1 ii = i - km1 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 ) ! !******************************************************************************* ! !! 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 ! untiary 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. ! ! Parameters: ! ! on entry ! ! r complex(ldr,p), where ldr >= p. ! r contains 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. ! ! Input, integer P, the order of the matrix. ! ! x complex(p). ! x contains the row to be added to r. x is ! not altered by cchud. ! ! z complex(ldz,nz), where ldz >= p. ! z is an array containing nz p-vectors to ! be updated with r. ! ! Input, integer LDZ, the leading dimension of Z. ! ! nz integer. ! nz is the number of vectors to be updated ! nz may be zero, in which case z, y, and rho ! are not referenced. ! ! y complex(nz). ! y contains the scalars for updating the vectors ! z. y is not altered by cchud. ! ! rho real(nz). ! rho contains the norms of the residual ! vectors that are to be updated. if rho(j) ! is negative, it is left unaltered. ! ! on return ! ! rc ! rho contain the updated quantities. ! z ! ! c real(p). ! c contains the cosines of the transforming ! rotations. ! ! s complex(p). ! s contains 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 integer jm1 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. ! jm1 = j - 1 do i = 1, jm1 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 = cabs(zeta) if ( azeta /= 0.0e0 .and. rho(j) >= 0.0e0 ) 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 ) ! !******************************************************************************* ! !! CGBCO factors a complex band matrix by gaussian ! elimination and estimates the condition of the matrix. ! ! 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 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, lda >= 5 and abd should contain ! ! * * * + + + , * = not used ! * * 13 24 35 46 , + = used for pivoting ! * 12 23 34 45 56 ! 11 22 33 44 55 66 ! 21 32 43 54 65 * ! ! Parameters: ! ! on entry ! ! abd complex(lda, n) ! 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 . ! see the comments below for details. ! ! 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. ! ! ml integer ! number of diagonals below the main diagonal. ! 0 <= ml < n . ! ! mu integer ! number of diagonals above the main diagonal. ! 0 <= mu < n . ! more efficient if ml <= mu . ! ! on return ! ! abd 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. ! ! ipvt integer(n) ! an integer vector of pivot indices. ! ! rcond real ! 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. ! ! z complex(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 kb 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.0e0 l = ml + 1 is = l + mu do j = 1, n anorm = max (anorm,scasum(l,abd(is,j),1)) if (is > ml + 1) is = is - 1 if (j <= mu) l = l + 1 if (j >= n - ml) l = l - 1 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 = (1.0e0,0.0e0) do j = 1, n z(j) = (0.0e0,0.0e0) end do m = ml + mu + 1 ju = 0 do k = 1, n if (cabs1(z(k)) /= 0.0e0) ek = csign1(ek,-z(k)) if ( cabs1(ek-z(k)) > cabs1(abd(m,k)) ) then s = cabs1(abd(m,k))/cabs1(ek-z(k)) call csscal(n,s,z,1) ek = cmplx(s,0.0e0)*ek end if wk = ek - z(k) wkm = -ek - z(k) s = cabs1(wk) sm = cabs1(wkm) if ( cabs1(abd(m,k)) /= 0.0e0 ) then wk = wk/conjg(abd(m,k)) wkm = wkm/conjg(abd(m,k)) else wk = (1.0e0,0.0e0) wkm = (1.0e0,0.0e0) 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.0e0/scasum(n,z,1) call csscal(n,s,z,1) ! ! Solve hermitian(L) * Y = W. ! do kb = 1, n k = n + 1 - kb lm = min(ml,n-k) if (k < n) z(k) = z(k) + cdotc(lm,abd(m+1,k),1,z(k+1),1) if ( cabs1(z(k)) > 1.0e0 ) then s = 1.0e0/cabs1(z(k)) call csscal(n,s,z,1) end if l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = 1.0e0 ! ! 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) call caxpy(lm,t,abd(m+1,k),1,z(k+1),1) if ( cabs1(z(k)) > 1.0e0 ) then s = 1.0e0/cabs1(z(k)) call csscal(n,s,z,1) ynorm = s*ynorm end if end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = s*ynorm ! ! Solve U * Z = W. ! do kb = 1, n k = n + 1 - kb if ( cabs1(z(k)) > cabs1(abd(m,k)) ) then s = cabs1(abd(m,k))/cabs1(z(k)) call csscal(n,s,z,1) ynorm = s*ynorm end if if (cabs1(abd(m,k)) /= 0.0e0) z(k) = z(k)/abd(m,k) if (cabs1(abd(m,k)) == 0.0e0) z(k) = (1.0e0,0.0e0) lm = min(k,m) - 1 la = m - lm lz = k - lm t = -z(k) call caxpy(lm,t,abd(la,k),1,z(lz),1) end do ! ! Make ZNORM = 1. ! s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = s*ynorm if ( anorm /= 0.0e0 ) then rcond = ynorm/anorm else rcond = 0.0e0 end if return end subroutine cgbdi ( abd, lda, n, ml, mu, ipvt, det ) ! !******************************************************************************* ! !! CGBDI computes the determinant of a band matrix factored by CGBCO or CGBFA. ! ! ! Discussion: ! ! If the inverse is needed, use CGBSL N times. ! ! Parameters: ! ! on entry ! ! abd complex(lda, n) ! the output from cgbco or cgbfa. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! ml integer ! number of diagonals below the main diagonal. ! ! mu integer ! number of diagonals above the main diagonal. ! ! ipvt integer(n) ! the pivot vector from cgbco or cgbfa. ! ! on return ! ! det complex(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 . ! 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 real, parameter :: ten = 10.0e0 ! m = ml + mu + 1 det(1) = (1.0e0,0.0e0) det(2) = (0.0e0,0.0e0) do i = 1, n if ( ipvt(i) /= i ) then det(1) = -det(1) end if det(1) = abd(m,i)*det(1) if ( cabs1(det(1)) == 0.0e0 ) then exit end if do while ( cabs1(det(1)) < 1.0e0 ) det(1) = cmplx(ten,0.0e0)*det(1) det(2) = det(2) - (1.0e0,0.0e0) end do do while ( cabs1(det(1)) >= ten ) det(1) = det(1)/cmplx(ten,0.0e0) det(2) = det(2) + (1.0e0,0.0e0) end do end do return end subroutine cgbfa ( abd, lda, n, ml, mu, ipvt, info ) ! !******************************************************************************* ! !! CGBFA factors a complex band matrix by elimination. ! ! 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. ! ! Parameters: ! ! on entry ! ! abd complex(lda, n) ! 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 . ! see the comments below for details. ! ! 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. ! ! ml integer ! number of diagonals below the main diagonal. ! 0 <= ml < n . ! ! mu integer ! number of diagonals above the main diagonal. ! 0 <= mu < n . ! more efficient if ml <= mu . ! on return ! ! abd 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. ! ! ipvt integer(n) ! an integer vector of pivot indices. ! ! info integer ! = 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) = (0.0e0,0.0e0) 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) = (0.0e0,0.0e0) end if ! ! Find L = pivot index. ! lm = min(ml,n-k) l = icamax(lm+1,abd(m,k),1) + m - 1 ipvt(k) = l + k - m ! ! Zero pivot implies this column already triangularized. ! if (cabs1(abd(l,k)) == 0.0e0) then info = k cycle end if ! ! Interchange if necessary. ! if ( l /= m ) then call c_swap ( abd(l,k), abd(m,k) ) end if ! ! Compute multipliers. ! t = -(1.0e0,0.0e0)/abd(m,k) call cscal(lm,t,abd(m+1,k),1) ! ! 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 call caxpy(lm,t,abd(m+1,k),1,abd(mm+1,j),1) end do end do ipvt(n) = n if ( cabs1 ( abd(m,n) ) == 0.0e0 ) then info = n end if return end subroutine cgbsl ( abd, lda, n, ml, mu, ipvt, b, job ) ! !******************************************************************************* ! !! 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 RCOND > 0.0 ! 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 ! ! Parameters: ! ! on entry ! ! abd complex(lda, n) ! the output from cgbco or cgbfa. ! ! Input, integer LDA, the leading dimension of ABD. ! ! Input, integer N, the order of the matrix. ! ! ml integer ! number of diagonals below the main diagonal. ! ! mu integer ! number of diagonals above the main diagonal. ! ! ipvt integer(n) ! the pivot vector from cgbco or cgbfa. ! ! Input/output, complex B(N). On input, the right hand side. ! On output, the solution. ! ! job integer ! = 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 kb 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 call caxpy(lm,t,abd(m+1,k),1,b(k+1),1) end do end if ! ! Now solve U * X = Y. ! do kb = 1, n k = n + 1 - kb b(k) = b(k)/abd(m,k) lm = min(k,m) - 1 la = m - lm lb = k - lm t = -b(k) call caxpy(lm,t,abd(la,k),1,b(lb),1) 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,k),1,b(lb),1) b(k) = (b(k) - t)/conjg(abd(m,k)) end do ! ! Now solve hermitian ( L ) * X = Y. ! if (ml /= 0) then do kb = 1, n-1 k = n - kb lm = min(ml,n-k) b(k) = b(k) + cdotc(lm,abd(m+1,k),1,b(k+1),1) l = ipvt(k) if ( l /= k ) then call c_swap ( b(l), b(k) ) end if end do end if end if return end subroutine cgeco ( a, lda, n, ipvt, rcond, z ) ! !******************************************************************************* ! !! CGECO factors a complex matrix by gaussian elimination ! and estimates the condition of the matrix. ! ! 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. ! ! Parameters: ! ! on entry ! ! a complex(lda, n) ! the matrix to be factored. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! on return ! ! a 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. ! ! ipvt integer(n) ! an integer vector of pivot indices. ! ! rcond real ! 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. ! ! z complex(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 kb 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.0e0 do j = 1, n anorm = max (anorm,scasum(n,a(1,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 = (1.0e0,0.0e0) z(1:n) = (0.0e0,0.0e0) do k = 1, n if (cabs1(z(k)) /= 0.0e0) ek = csign1(ek,-z(k)) if ( cabs1(ek-z(k)) > cabs1(a(k,k)) ) then s = cabs1(a(k,k))/cabs1(ek-z(k)) call csscal(n,s,z,1) ek = cmplx(s,0.0e0)*ek end if wk = ek - z(k) wkm = -ek - z(k) s = cabs1(wk) sm = cabs1(wkm) if ( cabs1(a(k,k)) /= 0.0e0 ) then wk = wk/conjg(a(k,k)) wkm = wkm/conjg(a(k,k)) else wk = (1.0e0,0.0e0) wkm = (1.0e0,0.0e0) 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 do j = k+1, n z(j) = z(j) + t*conjg(a(k,j)) end do end if z(k) = wk end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ! ! Solve hermitian(l)*y = W. ! do kb = 1, n k = n + 1 - kb if (k < n) z(k) = z(k) + cdotc(n-k,a(k+1,k),1,z(k+1),1) if ( cabs1(z(k)) > 1.0e0 ) then s = 1.0e0/cabs1(z(k)) call csscal(n,s,z,1) end if l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = 1.0e0 ! ! Solve l*v = Y. ! do k = 1, n l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t if (k < n) call caxpy(n-k,t,a(k+1,k),1,z(k+1),1) if ( cabs1(z(k)) > 1.0e0 ) then s = 1.0e0/cabs1(z(k)) call csscal(n,s,z,1) ynorm = s*ynorm end if end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = s*ynorm ! ! Solve u*z = V. ! do kb = 1, n k = n + 1 - kb if ( cabs1(z(k)) > cabs1(a(k,k)) ) then s = cabs1(a(k,k))/cabs1(z(k)) call csscal(n,s,z,1) ynorm = s*ynorm end if if (cabs1(a(k,k)) /= 0.0e0) z(k) = z(k)/a(k,k) if (cabs1(a(k,k)) == 0.0e0) z(k) = (1.0e0,0.0e0) t = -z(k) call caxpy(k-1,t,a(1,k),1,z(1),1) end do ! ! Make ZNORM = 1. ! s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = s*ynorm if ( anorm /= 0.0e0 ) then rcond = ynorm / anorm else rcond = 0.0e0 end if return end subroutine cgedi ( a, lda, n, ipvt, det, work, job ) ! !******************************************************************************* ! !! 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 rcond > 0.0 or cgefa has set ! info == 0 . ! ! Parameters: ! ! on entry ! ! a complex(lda, n) ! the output from cgeco or cgefa. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! ipvt integer(n) ! the pivot vector from cgeco or cgefa. ! ! work complex(n) ! work vector. contents destroyed. ! ! job integer ! = 11 both determinant and inverse. ! = 01 inverse only. ! = 10 determinant only. ! ! on return ! ! a inverse of original matrix if requested. ! otherwise unchanged. ! ! det complex(2) ! determinant of 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 . ! 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 kb integer kp1 integer l complex t real, parameter :: ten = 10.0e0 complex work(n) ! ! Compute the determinant. ! if ( job / 10 /= 0 ) then det(1) = (1.0e0,0.0e0) det(2) = (0.0e0,0.0e0) do i = 1, n if (ipvt(i) /= i) det(1) = -det(1) det(1) = a(i,i)*det(1) if ( cabs1(det(1)) == 0.0e0 ) then exit end if do while ( cabs1(det(1)) < 1.0e0 ) det(1) = cmplx(ten,0.0e0)*det(1) det(2) = det(2) - (1.0e0,0.0e0) end do do while ( cabs1(det(1)) >= ten ) det(1) = det(1)/cmplx(ten,0.0e0) det(2) = det(2) + (1.0e0,0.0e0) end do end do end if ! ! Compute inverse(U). ! if ( mod(job,10) /= 0 ) then do k = 1, n a(k,k) = (1.0e0,0.0e0)/a(k,k) t = -a(k,k) call cscal(k-1,t,a(1,k),1) kp1 = k + 1 do j = k+1, n t = a(k,j) a(k,j) = (0.0e0,0.0e0) call caxpy(k,t,a(1,k),1,a(1,j),1) end do end do ! ! Form inverse(u) * inverse(l). ! do kb = 1, n-1 k = n - kb kp1 = k + 1 do i = k+1, n work(i) = a(i,k) a(i,k) = (0.0e0,0.0e0) end do do j = k+1, n t = work(j) call caxpy(n,t,a(1,j),1,a(1,k),1) end do l = ipvt(k) if (l /= k) call cswap(n,a(1,k),1,a(1,l),1) end do end if return end subroutine cgefa ( a, lda, n, ipvt, info ) ! !******************************************************************************* ! !! CGEFA factors a complex matrix by Gaussian elimination. ! ! ! Parameters: ! ! on entry ! ! a complex(lda, n) ! the matrix to be factored. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! on return ! ! a 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. ! ! ipvt integer(n) ! an integer vector of pivot indices. ! ! info integer ! = 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 kp1 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,k), 1 ) + k - 1 ipvt(k) = l ! ! Zero pivot implies this column already triangularized. ! if ( cabs1 ( a(l,k) ) == 0.0e0 ) then info = k cycle end if ! ! Interchange if necessary. ! if ( l /= k ) then call c_swap ( a(l,k), a(k,k) ) end if ! ! Compute multipliers ! t = -( 1.0e0, 0.0e0 ) / a(k,k) call cscal ( n-k, t, a(k+1,k), 1 ) ! ! 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 call caxpy ( n-k, t, a(k+1,k), 1, a(k+1,j), 1 ) end do end do ipvt(n) = n if ( cabs1 ( a(n,n) ) == 0.0e0 ) then info = n end if return end subroutine cgesl ( a, lda, n, ipvt, b, job ) ! !******************************************************************************* ! !! 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 rcond > 0.0 ! 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 ! ! 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. ! ! ipvt 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. ! ! job 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 kb 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 call caxpy ( n-k, t, a(k+1,k), 1, b(k+1), 1 ) end do ! ! Now solve u*x = y. ! do kb = 1, n k = n + 1 - kb b(k) = b(k) / a(k,k) t = -b(k) call caxpy ( k-1, t, a(1,k), 1, b(1), 1 ) 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, b(1), 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,k), 1, b(k+1), 1 ) 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 ) ! !******************************************************************************* ! !! CGTSL solves a complex general tridiagonal system. ! ! ! Parameters: ! ! on entry ! ! Input, integer N, the order of the matrix. ! ! c complex(n) ! is the subdiagonal of the tridiagonal matrix. ! c(2) through c(n) should contain the subdiagonal. ! on output c is destroyed. ! ! d complex(n) ! is the diagonal of the tridiagonal matrix. ! on output d is destroyed. ! ! e complex(n) ! is the superdiagonal of the tridiagonal matrix. ! e(1) through e(n-1) should contain the superdiagonal. ! on output e is destroyed. ! ! Input/output, complex B(N). On input, the right hand side. ! On output, the solution. ! ! on return ! ! info integer ! = 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 ( n-1 >= 1 ) then d(1) = e(1) e(1) = ( 0.0e0, 0.0e0 ) e(n) = ( 0.0e0, 0.0e0 ) do k = 1, n-1 if ( cabs1 ( c(k+1) ) >= cabs1 ( c(k) ) ) then call c_swap ( c(k+1), c(k) ) call c_swap ( d(k+1), d(k) ) call c_swap ( e(k+1), e(k) ) call c_swap ( b(k+1), b(k) ) end if if ( cabs1 ( c(k) ) == 0.0e0 ) 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) = ( 0.0e0, 0.0e0 ) b(k+1) = b(k+1) + t * b(k) end do end if if ( cabs1 ( c(n) ) == 0.0e0 ) then info = n return end if ! ! Back solve. ! b(n) = b(n) / c(n) if ( n > 1 ) 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 ) ! !******************************************************************************* ! !! CHICO factors a complex hermitian matrix by elimination with ! symmetric pivoting and estimates the condition of the matrix. ! ! 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. ! ! Parameters: ! ! on entry ! ! a complex(lda, n) ! the hermitian matrix to be factored. ! 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 ! ! a 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. ! ! ipvt integer(n) ! an integer vector of pivot indices. ! ! rcond real ! 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. ! ! z complex(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 jm1 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),1),0.0e0) jm1 = j - 1 do i = 1, jm1 z(i) = cmplx(real(z(i))+cabs1(a(i,j)),0.0e0) end do end do anorm = 0.0e0 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 = (1.0e0,0.0e0) z(1:n) = (0.0e0,0.0e0) k = n do while ( k > 0 ) ks = 1 if (ipvt(k) < 0) ks = 2 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.0e0) ek = csign1(ek,z(k)) z(k) = z(k) + ek call caxpy(k-ks,z(k),a(1,k),1,z(1),1) if ( ks /= 1 ) then if (cabs1(z(k-1)) /= 0.0e0) ek = csign1(ek,z(k-1)) z(k-1) = z(k-1) + ek call caxpy(k-ks,z(k-1),a(1,k-1),1,z(1),1) end if if ( ks /= 2 ) then if ( cabs1(z(k)) > cabs1(a(k,k)) ) then s = cabs1(a(k,k))/cabs1(z(k)) call csscal(n,s,z,1) ek = cmplx(s,0.0e0)*ek end if if (cabs1(a(k,k)) /= 0.0e0) z(k) = z(k)/a(k,k) if (cabs1(a(k,k)) == 0.0e0) z(k) = (1.0e0,0.0e0) 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.0e0 z(k) = (akm1*bk - bkm1)/denom z(k-1) = (ak*bkm1 - bk)/denom end if k = k - ks end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ! ! Solve hermitian(u)*y = w. ! k = 1 do while ( k <= n ) ks = 1 if (ipvt(k) < 0) ks = 2 if ( k /= 1 ) then z(k) = z(k) + cdotc(k-1,a(1,k),1,z(1),1) if (ks == 2) then z(k+1) = z(k+1) + cdotc(k-1,a(1,k+1),1,z(1),1) end if kp = abs(ipvt(k)) if ( kp /= k ) then call c_swap ( z(k), z(kp) ) end if end if k = k + ks end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = 1.0e0 ! ! Solve u*d*v = y. ! k = n do while ( k > 0 ) ks = 1 if (ipvt(k) < 0) ks = 2 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 call caxpy(k-ks,z(k),a(1,k),1,z(1),1) if (ks == 2) call caxpy(k-ks,z(k-1),a(1,k-1),1,z(1),1) end if if ( ks /= 2 ) then if ( cabs1(z(k)) > cabs1(a(k,k)) ) then s = cabs1(a(k,k))/cabs1(z(k)) call csscal(n,s,z,1) ynorm = s*ynorm end if if (cabs1(a(k,k)) /= 0.0e0) z(k) = z(k)/a(k,k) if (cabs1(a(k,k)) == 0.0e0) z(k) = (1.0e0,0.0e0) 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.0e0 z(k) = (akm1*bk - bkm1)/denom z(k-1) = (ak*bkm1 - bk)/denom end if k = k - ks end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = s*ynorm ! ! Solve hermitian(u)*z = v. ! k = 1 do while ( k <= n ) ks = 1 if (ipvt(k) < 0) ks = 2 if ( k /= 1 ) then z(k) = z(k) + cdotc(k-1,a(1,k),1,z(1),1) if (ks == 2) & z(k+1) = z(k+1) + cdotc(k-1,a(1,k+1),1,z(1),1) kp = abs(ipvt(k)) if ( kp /= k ) then call c_swap ( z(k), z(kp) ) end if end if k = k + ks end do ! ! Make ZNORM = 1. ! s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = s*ynorm if ( anorm /= 0.0e0 ) then rcond = ynorm/anorm else rcond = 0.0e0 end if return end subroutine chidi ( a, lda, n, ipvt, det, inert, work, job ) ! !******************************************************************************* ! !! 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 . ! ! Parameters: ! ! on entry ! ! a complex(lda,n) ! the output from chifa. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! ipvt integer(n) ! the pivot vector from chifa. ! ! work complex(n) ! work vector. contents destroyed. ! ! job 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. ! ! on return ! ! variables not requested by job are not used. ! ! a contains the upper triangle of the inverse of ! the original matrix. the strict lower triangle ! is never referenced. ! ! det real(2) ! determinant of original matrix. ! determinant = det(1) * 10.0**det(2) ! with 1.0 <= abs(det(1)) < 10.0 ! or det(1) = 0.0. ! ! inert integer(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. ! 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 real, parameter :: ten = 10.0e0 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) = 0 inert(2) = 0 inert(3) = 0 end if if ( .not. nodet ) then det(1) = 1.0e0 det(2) = 0.0e0 end if t = 0.0e0 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 = cabs(s) ! to avoid underflow/overflow troubles. ! take two passes through scaling. Use T for flag. ! if ( t == 0.0e0 ) then t = cabs(a(k,k+1)) d = (d/t)*real(a(k+1,k+1)) - t else d = t t = 0.0e0 end if end if if ( .not. noert ) then if ( d > 0.0e0 ) then inert(1) = inert(1) + 1 else if ( d < 0.0e0 ) then inert(2) = inert(2) + 1 else if ( d == 0.0e0 ) then inert(3) = inert(3) + 1 end if end if if ( .not. nodet ) then det(1) = d*det(1) if ( det(1) /= 0.0e0 ) then do while (abs(det(1)) < 1.0e0) det(1) = ten*det(1) det(2) = det(2) - 1.0e0 end do do while ( abs(det(1)) >= ten ) det(1) = det(1)/ten det(2) = det(2) + 1.0e0 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 ( ipvt(k) >= 0 ) then ! ! 1 by 1 ! a(k,k) = cmplx(1.0e0/real(a(k,k)),0.0e0) if ( km1 >= 1 ) then call ccopy(km1,a(1,k),1,work,1) do j = 1, km1 a(j,k) = cdotc(j,a(1,j),1,work,1) call caxpy(j-1,work(j),a(1,j),1,a(1,k),1) end do a(k,k) = a(k,k) & + cmplx(real(cdotc(km1,work,1,a(1,k),1)),0.0e0) end if kstep = 1 else ! ! 2 by 2 ! t = cabs(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.0e0) a(k,k) = cmplx(akp1/d,0.0e0) a(k+1,k+1) = cmplx(ak/d,0.0e0) a(k,k+1) = -akkp1/d if ( km1 >= 1 ) then call ccopy(km1,a(1,k+1),1,work,1) do j = 1, km1 a(j,k+1) = cdotc(j,a(1,j),1,work,1) call caxpy(j-1,work(j),a(1,j),1,a(1,k+1),1) end do a(k+1,k+1) = a(k+1,k+1) & + cmplx(real(cdotc(km1,work,1,a(1,k+1),1)),0.0e0) a(k,k+1) = a(k,k+1) + cdotc(km1,a(1,k),1,a(1,k+1),1) call ccopy(km1,a(1,k),1,work,1) do j = 1, km1 a(j,k) = cdotc(j,a(1,j),1,work,1) call caxpy(j-1,work(j),a(1,j),1,a(1,k),1) end do a(k,k) = a(k,k) & + cmplx(real(cdotc(km1,work,1,a(1,k),1)),0.0e0) end if kstep = 2 end if ! ! Swap ! ks = abs(ipvt(k)) if ( ks /= k ) then call cswap(ks,a(1,ks),1,a(1,k),1) do jb = ks, k j = k + ks - jb temp = conjg(a(j,k)) a(j,k) = conjg(a(ks,j)) a(ks,j) = temp end do if ( kstep /= 1 ) then call c_swap ( a(ks,k+1), a(k,k+1) ) end if end if k = k + kstep end do end if return end subroutine chifa ( a, lda, n, ipvt, info ) ! !******************************************************************************* ! !! 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. ! ! Parameters: ! ! on entry ! ! a complex(lda,n) ! the hermitian matrix to be factored. ! only the diagonal and upper triangle are used. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! on return ! ! a 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. ! ! ipvt integer(n) ! an integer vector of pivot indices. ! ! info integer ! = 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 imaxp1 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.0e0 + sqrt(17.0e0))/8.0e0 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.0e0) info = 1 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) colmax = cabs1(a(imax,k)) if ( absakk >= alpha*colmax) then kstep = 1 swap = .false. else ! ! Determine the largest off-diagonal element in row IMAX. ! rowmax = 0.0e0 imaxp1 = imax + 1 do j = imaxp1, k rowmax = max (rowmax,cabs1(a(imax,j))) end do if ( imax /= 1 ) then jmax = icamax(imax-1,a(1,imax),1) rowmax = max (rowmax,cabs1(a(jmax,imax))) end if if ( cabs1(a(imax,imax)) >= alpha*rowmax ) then kstep = 1 swap = .true. else if ( absakk >= alpha*colmax*(colmax/rowmax)) 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.0e0) 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),1,a(1,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) call caxpy(j,t,a(1,k),1,a(1,j),1) a(j,j) = cmplx(real(a(j,j)),0.0e0) a(j,k) = mulk end do ! ! Set the pivot array. ! ipvt(k) = k if (swap) ipvt(k) = imax else ! ! 2 x 2 pivot block. ! if ( swap ) then call cswap(imax,a(1,imax),1,a(1,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 ( km2 > 0 ) then ak = a(k,k)/a(k-1,k) akm1 = a(k-1,k-1)/conjg(a(k-1,k)) denom = 1.0e0 - ak*akm1 do jj = 1, km2 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) call caxpy(j,t,a(1,k),1,a(1,j),1) t = conjg(mulkm1) call caxpy(j,t,a(1,k-1),1,a(1,j),1) a(j,k) = mulk a(j,k-1) = mulkm1 a(j,j) = cmplx(real(a(j,j)),0.0e0) end do end if ! ! Set the pivot array. ! ipvt(k) = 1 - k if (swap) ipvt(k) = -imax ipvt(k-1) = ipvt(k) end if k = k - kstep end do return end subroutine chisl ( a, lda, n, ipvt, b ) ! !******************************************************************************* ! !! 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 ! ! Parameters: ! ! on entry ! ! a complex(lda,n) ! the output from chifa. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! ipvt integer(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 temp ! ! Loop backward applying the transformations and d inverse to b. ! k = n do while ( k > 0 ) ! ! 1 x 1 pivot block. ! if ( ipvt(k) >= 0 ) then if ( k /= 1 ) then kp = ipvt(k) if ( kp /= k ) then call c_swap ( b(k), b(kp) ) end if call caxpy(k-1,b(k),a(1,k),1,b(1),1) 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 call c_swap ( b(k-1), b(kp) ) end if call caxpy(k-2,b(k),a(1,k),1,b(1),1) call caxpy(k-2,b(k-1),a(1,k-1),1,b(1),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.0e0 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 (ipvt(k) >= 0) then if ( k /= 1 ) then b(k) = b(k) + cdotc(k-1,a(1,k),1,b(1),1) kp = ipvt(k) if ( kp /= k ) then temp = b(k) b(k) = b(kp) b(kp) = temp 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,b(1),1) b(k+1) = b(k+1) + cdotc(k-1,a(1,k+1),1,b(1),1) kp = abs(ipvt(k)) if ( kp /= k ) then call c_swap ( b(k), b(kp) ) end if end if k = k + 2 end if end do return end subroutine chpco ( ap, n, ipvt, rcond, z ) ! !******************************************************************************* ! !! CHPCO factors a complex hermitian matrix stored in packed ! form by elimination with symmetric pivoting and estimates ! the condition of the matrix. ! ! 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 ! ! Parameters: ! ! on entry ! ! ap complex (n*(n+1)/2) ! 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 . ! see comments below for details. ! ! Input, integer N, the order of the matrix. ! ! output ! ! ap 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. ! ! ipvt integer(n) ! an integer vector of pivot indices. ! ! rcond real ! 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. ! ! z complex(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 jm1 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),1),0.0e0) ij = j1 j1 = j1 + j jm1 = j - 1 do i = 1, jm1 z(i) = cmplx(real(z(i))+cabs1(ap(ij)),0.0e0) ij = ij + 1 end do end do anorm = 0.0e0 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 = (1.0e0,0.0e0) z(1:n) = (0.0e0,0.0e0) k = n ik = (n*(n - 1))/2 do while ( k > 0 ) kk = ik + k ikm1 = ik - (k - 1) ks = 1 if (ipvt(k) < 0) ks = 2 kp = abs(ipvt(k)) kps = k + 1 - ks if ( kp /= kps ) then call c_swap ( z(kps), z(kp) ) end if if (cabs1(z(k)) /= 0.0e0) ek = csign1(ek,z(k)) z(k) = z(k) + ek call caxpy(k-ks,z(k),ap(ik+1),1,z(1),1) if ( ks /= 1 ) then if (cabs1(z(k-1)) /= 0.0e0) ek = csign1(ek,z(k-1)) z(k-1) = z(k-1) + ek call caxpy(k-ks,z(k-1),ap(ikm1+1),1,z(1),1) end if if ( ks /= 2 ) then if ( cabs1(z(k)) > cabs1(ap(kk)) ) then s = cabs1(ap(kk))/cabs1(z(k)) call csscal(n,s,z,1) ek = cmplx(s,0.0e0)*ek end if if (cabs1(ap(kk)) /= 0.0e0) z(k) = z(k)/ap(kk) if (cabs1(ap(kk)) == 0.0e0) z(k) = (1.0e0,0.0e0) 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.0e0 z(k) = (akm1*bk - bkm1)/denom z(k-1) = (ak*bkm1 - bk)/denom end if k = k - ks ik = ik - k if (ks == 2) ik = ik - (k + 1) end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ! ! Solve hermitian(u)*y = w. ! k = 1 ik = 0 do while ( k <= n ) ks = 1 if (ipvt(k) < 0) ks = 2 if ( k /= 1 ) then z(k) = z(k) + cdotc(k-1,ap(ik+1),1,z(1),1) ikp1 = ik + k if (ks == 2) then z(k+1) = z(k+1) + cdotc(k-1,ap(ikp1+1),1,z(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) ik = ik + (k + 1) k = k + ks end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = 1.0e0 ! ! Solve u*d*v = y. ! k = n ik = n*(n - 1)/2 do while ( k > 0 ) kk = ik + k ikm1 = ik - (k - 1) ks = 1 if (ipvt(k) < 0) ks = 2 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 call caxpy(k-ks,z(k),ap(ik+1),1,z(1),1) if (ks == 2) call caxpy(k-ks,z(k-1),ap(ikm1+1),1,z(1),1) end if if ( ks /= 2 ) then if ( cabs1(z(k)) > cabs1(ap(kk))) then s = cabs1(ap(kk))/cabs1(z(k)) call csscal(n,s,z,1) ynorm = s*ynorm end if if (cabs1(ap(kk)) /= 0.0e0) z(k) = z(k)/ap(kk) if (cabs1(ap(kk)) == 0.0e0) z(k) = (1.0e0,0.0e0) 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.0e0 z(k) = (akm1*bk - bkm1)/denom z(k-1) = (ak*bkm1 - bk)/denom end if k = k - ks ik = ik - k if (ks == 2) ik = ik - (k + 1) end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = s*ynorm ! ! Solve hermitian(u)*z = v. ! k = 1 ik = 0 do while ( k <= n ) ks = 1 if (ipvt(k) < 0) ks = 2 if (k /= 1) then z(k) = z(k) + cdotc(k-1,ap(ik+1),1,z(1),1) ikp1 = ik + k if (ks == 2) then z(k+1) = z(k+1) + cdotc(k-1,ap(ikp1+1),1,z(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) ik = ik + (k + 1) k = k + ks end do ! ! Make ZNORM = 1. ! s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = s*ynorm if ( anorm /= 0.0e0 ) then rcond = ynorm/anorm else rcond = 0.0e0 end if return end subroutine chpdi ( ap, n, ipvt, det, inert, work, job ) ! !******************************************************************************* ! !! CHPDI computes the determinant, inertia and inverse ! of a complex hermitian matrix using the factors from chpfa, ! where the matrix is stored in packed form. ! ! error condition ! ! a division by zero will occur if the inverse is requested ! and chpco has set rcond == 0.0 ! or chpfa has set info /= 0 . ! ! Parameters: ! ! on entry ! ! ap complex (n*(n+1)/2) ! the output from chpfa. ! ! Input, integer N, the order of the matrix. ! ! ipvt integer(n) ! the pivot vector from chpfa. ! ! work complex(n) ! work vector. contents ignored. ! ! job 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. ! ! on return ! ! variables not requested by job are not used. ! ! 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. ! ! det real(2) ! determinant of original matrix. ! determinant = det(1) * 10.0**det(2) ! with 1.0 <= abs(det(1)) < 10.0 ! or det(1) = 0.0. ! ! inert integer(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. ! 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 real, parameter :: ten = 10.0E+00 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) = 0 inert(2) = 0 inert(3) = 0 end if if ( .not. nodet ) then det(1) = 1.0e0 det(2) = 0.0e0 end if t = 0.0e0 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 = cabs(s) ! to avoid underflow/overflow troubles. ! take two passes through scaling. use t for flag. ! if ( t == 0.0e0 ) then ikp1 = ik + k kkp1 = ikp1 + k t = cabs(ap(kkp1)) d = (d/t)*real(ap(kkp1+1)) - t else d = t t = 0.0e0 end if end if if ( .not. noert ) then if (d > 0.0e0) inert(1) = inert(1) + 1 if (d < 0.0e0) inert(2) = inert(2) + 1 if (d == 0.0e0) inert(3) = inert(3) + 1 end if if ( .not. nodet ) then det(1) = d*det(1) if (det(1) /= 0.0e0 ) then do while (abs(det(1)) < 1.0e0) det(1) = ten*det(1) det(2) = det(2) - 1.0e0 end do do while (abs(det(1)) >= ten) det(1) = det(1)/ten det(2) = det(2) + 1.0e0 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 ( ipvt(k) >= 0 ) then ap(kk) = cmplx(1.0e0/real(ap(kk)),0.0e0) if ( km1 >= 1 ) then call ccopy(km1,ap(ik+1),1,work,1) ij = 0 do j = 1, km1 jk = ik + j ap(jk) = cdotc(j,ap(ij+1),1,work,1) call caxpy(j-1,work(j),ap(ij+1),1,ap(ik+1),1) ij = ij + j end do ap(kk) = ap(kk) & + cmplx(real(cdotc(km1,work,1,ap(ik+1),1)), 0.0e0) end if kstep = 1 ! ! 2 by 2 ! else t = cabs(ap(kkp1)) ak = real(ap(kk))/t akp1 = real(ap(kkp1+1))/t akkp1 = ap(kkp1)/t d = t*(ak*akp1 - 1.0e0) ap(kk) = cmplx(akp1/d,0.0e0) ap(kkp1+1) = cmplx(ak/d,0.0e0) ap(kkp1) = -akkp1/d if ( km1 >= 1 ) then call ccopy(km1,ap(ikp1+1),1,work,1) ij = 0 do j = 1, km1 jkp1 = ikp1 + j ap(jkp1) = cdotc(j,ap(ij+1),1,work,1) call caxpy(j-1,work(j),ap(ij+1),1,ap(ikp1+1),1) ij = ij + j end do ap(kkp1+1) = ap(kkp1+1) + cmplx(real(cdotc(km1,work,1, & ap(ikp1+1),1)),0.0e0) ap(kkp1) = ap(kkp1) + cdotc(km1,ap(ik+1),1,ap(ikp1+1),1) call ccopy(km1,ap(ik+1),1,work,1) ij = 0 do j = 1, km1 jk = ik + j ap(jk) = cdotc(j,ap(ij+1),1,work,1) call caxpy(j-1,work(j),ap(ij+1),1,ap(ik+1),1) ij = ij + j end do ap(kk) = ap(kk) & + cmplx(real(cdotc(km1,work,1,ap(ik+1),1)), 0.0e0) 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),1,ap(ik+1),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) ik = ik + k + 1 k = k + kstep end do end if return end subroutine chpfa ( ap, n, ipvt, info ) ! !******************************************************************************* ! !! CHPFA factors a complex hermitian matrix stored in ! packed form by elimination with symmetric pivoting. ! ! 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 ! ! Parameters: ! ! on entry ! ! ap complex (n*(n+1)/2) ! 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 . ! see comments below for details. ! ! Input, integer N, the order of the matrix. ! ! output ! ! ap 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. ! ! ipvt integer(n) ! an integer vector of pivot indices. ! ! info integer ! = 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 imaxp1 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.0e0 + sqrt(17.0e0))/8.0e0 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.0e0) info = 1 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),1) imk = ik + imax colmax = cabs1(ap(imk)) if ( absakk >= alpha*colmax) then kstep = 1 swap = .false. ! ! Determine the largest off-diagonal element in row IMAX. ! else rowmax = 0.0e0 imaxp1 = imax + 1 im = imax*(imax - 1)/2 imj = im + 2*imax do j = imaxp1, k rowmax = max (rowmax,cabs1(ap(imj))) imj = imj + j end do if ( imax /= 1 ) then jmax = icamax(imax-1,ap(im+1),1) jmim = jmax + im rowmax = max (rowmax,cabs1(ap(jmim))) end if imim = imax + im if ( cabs1(ap(imim)) >= alpha*rowmax) then kstep = 1 swap = .true. else if ( absakk >= alpha*colmax*(colmax/rowmax) ) 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.0e0) then ipvt(k) = k info = k ik = ik - (k - 1) if (kstep == 2) ik = ik - (k - 2) k = k - kstep cycle end if if ( kstep /= 2 ) then ! ! 1 x 1 pivot block. ! if ( swap ) then call cswap(imax,ap(im+1),1,ap(ik+1),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) call caxpy(j,t,ap(ik+1),1,ap(ij+1),1) ijj = ij + j ap(ijj) = cmplx(real(ap(ijj)),0.0e0) ap(jk) = mulk ij = ij - (j - 1) end do ! ! Set the pivot array. ! ipvt(k) = k if (swap) ipvt(k) = imax ! ! 2 x 2 pivot block. ! else km1k = ik + k - 1 ikm1 = ik - (k - 1) if ( swap ) then call cswap(imax,ap(im+1),1,ap(ikm1+1),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.0e0 - 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) call caxpy(j,t,ap(ik+1),1,ap(ij+1),1) t = conjg(mulkm1) call caxpy(j,t,ap(ikm1+1),1,ap(ij+1),1) ap(jk) = mulk ap(jkm1) = mulkm1 ijj = ij + j ap(ijj) = cmplx(real(ap(ijj)),0.0e0) ij = ij - (j - 1) end do end if ! ! Set the pivot array. ! ipvt(k) = 1 - k if (swap) ipvt(k) = -imax ipvt(k-1) = ipvt(k) end if ik = ik - (k - 1) if (kstep == 2) ik = ik - (k - 2) k = k - kstep end do return end subroutine chpsl ( ap, n, ipvt, b ) ! !******************************************************************************* ! !! 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 ! ! Parameters: ! ! on entry ! ! ap complex(n*(n+1)/2) ! the output from chpfa. ! ! Input, integer N, the order of the matrix. ! ! ipvt integer(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 temp ! ! Loop backward applying the transformations and inverse ( D ) to B. ! k = n ik = (n*(n - 1))/2 do while ( k > 0 ) kk = ik + k ! ! 1 x 1 pivot block. ! if (ipvt(k) >= 0 ) then if ( k /= 1 ) then kp = ipvt(k) if ( kp /= k ) then call c_swap ( b(k), b(kp) ) end if call caxpy(k-1,b(k),ap(ik+1),1,b(1),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 call c_swap ( b(k-1), b(kp) ) end if call caxpy(k-2,b(k),ap(ik+1),1,b(1),1) call caxpy(k-2,b(k-1),ap(ikm1+1),1,b(1),1) 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.0e0 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 ( ipvt(k) >= 0) then if ( k /= 1) then b(k) = b(k) + cdotc(k-1,ap(ik+1),1,b(1),1) kp = ipvt(k) if ( kp /= k ) then temp = b(k) b(k) = b(kp) b(kp) = temp 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),1,b(1),1) ikp1 = ik + k b(k+1) = b(k+1) + cdotc(k-1,ap(ikp1+1),1,b(1),1) kp = abs(ipvt(k)) if ( kp /= k ) then call c_swap ( b(k), b(kp) ) 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 ) ! !******************************************************************************* ! !! CPBCO factors a complex hermitian positive definite matrix ! stored in band form and estimates the condition 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 ! ! Parameters: ! ! on entry ! ! abd complex(lda, n) ! 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 . see the comments below for details. ! ! Input, integer LDA, the leading dimension of ABD. ! LDA must be at least M+1. ! ! Input, integer N, the order of the matrix. ! ! m integer ! the number of diagonals above the main diagonal. ! 0 <= m < n . ! ! on return ! ! abd an upper triangular matrix r , stored in band ! form, so that a = hermitian(r)*r . ! if info /= 0 , the factorization is not complete. ! ! rcond real ! 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. if info /= 0 , rcond is unchanged. ! ! z complex(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. ! ! info integer ! = 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 kb 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) ! ! 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,j),1),0.0e0) k = j - l do i = mu, m k = k + 1 z(k) = cmplx(real(z(k))+cabs1(abd(i,j)),0.0e0) end do end do anorm = 0.0e0 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 = (1.0e0,0.0e0) z(1:n) = (0.0e0,0.0e0) do k = 1, n if (cabs1(z(k)) /= 0.0e0) ek = csign1(ek,-z(k)) if ( cabs1(ek-z(k)) > real(abd(m+1,k)) ) then s = real(abd(m+1,k))/cabs1(ek-z(k)) call csscal(n,s,z,1) ek = cmplx(s,0.0e0)*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) kp1 = k + 1 j2 = min(k+m,n) i = m + 1 if ( kp1 <= j2 ) then do j = kp1, 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 = kp1, 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.0e0/scasum(n,z,1) call csscal(n,s,z,1) ! ! Solve r*y = w ! do kb = 1, n k = n + 1 - kb if ( cabs1(z(k)) > real(abd(m+1,k)) ) then s = real(abd(m+1,k))/cabs1(z(k)) call csscal(n,s,z,1) 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) call caxpy(lm,t,abd(la,k),1,z(lb),1) end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = 1.0e0 ! ! 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,k),1,z(lb),1) if ( cabs1(z(k)) > real(abd(m+1,k)) ) then s = real(abd(m+1,k))/cabs1(z(k)) call csscal(n,s,z,1) ynorm = s*ynorm end if z(k) = z(k)/abd(m+1,k) end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = s*ynorm ! ! Solve r*z = w. ! do kb = 1, n k = n + 1 - kb if ( cabs1(z(k)) > real(abd(m+1,k)) ) then s = real(abd(m+1,k))/cabs1(z(k)) call csscal(n,s,z,1) 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) call caxpy(lm,t,abd(la,k),1,z(lb),1) end do ! ! Make ZNORM = 1. ! s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = s*ynorm if ( anorm /= 0.0e0 ) then rcond = ynorm/anorm else rcond = 0.0e0 end if return end subroutine cpbdi ( abd, lda, n, m, det ) ! !******************************************************************************* ! !! 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. ! ! Parameters: ! ! on entry ! ! 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. ! ! m integer ! the number of diagonals above the main diagonal. ! ! on return ! ! det real(2) ! determinant of 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 real, parameter :: ten = 10.0E+00 ! det(1) = 1.0e0 det(2) = 0.0e0 do i = 1, n det(1) = real(abd(m+1,i))**2 * det(1) if ( det(1) == 0.0e0 ) then exit end if do while ( det(1) < 1.0e0 ) det(1) = ten * det(1) det(2) = det(2) - 1.0e0 end do do while ( det(1) >= ten ) det(1) = det(1) / ten det(2) = det(2) + 1.0e0 end do end do return end subroutine cpbfa ( abd, lda, n, m, info ) ! !******************************************************************************* ! !! 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 ! ! Parameters: ! ! on entry ! ! abd complex(lda, n) ! 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 . see the comments below for details. ! ! Input, integer LDA, the leading dimension of ABD. ! LDA must be at least M+1. ! ! Input, integer N, the order of the matrix. ! ! m integer ! the number of diagonals above the main diagonal. ! 0 <= m < n . ! ! on return ! ! abd an upper triangular matrix r , stored in band ! form, so that a = hermitian(r)*r . ! ! info integer ! = 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 info = j s = 0.0e0 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,jk),1,abd(mu,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.0e0 .or. aimag(abd(m+1,j)) /= 0.0e0) then info = j exit end if abd(m+1,j) = cmplx(sqrt(s),0.0e0) end do return end subroutine cpbsl ( abd, lda, n, m, b ) ! !******************************************************************************* ! !! CPBSL solves the complex hermitian positive definite band ! system a*x = b ! using the factors computed by cpbco or cpbfa. ! ! error condition ! ! 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) go to ... ! do j = 1, p ! call cpbsl(abd,lda,n,c(1,j)) ! end do ! ! Parameters: ! ! on entry ! ! abd complex(lda, n) ! the output from cpbco or cpbfa. ! ! Input, integer LDA, the leading dimension of ABD. ! ! Input, integer N, the order of the matrix. ! ! m integer ! 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 kb 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,k),1,b(lb),1) b(k) = (b(k) - t)/abd(m+1,k) end do ! ! Solve r*x = y ! do kb = 1, n k = n + 1 - kb lm = min(k-1,m) la = m + 1 - lm lb = k - lm b(k) = b(k)/abd(m+1,k) t = -b(k) call caxpy(lm,t,abd(la,k),1,b(lb),1) end do return end subroutine cpoco ( a, lda, n, rcond, z, info ) ! !******************************************************************************* ! !! CPOCO factors a complex hermitian positive definite matrix ! and 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. ! ! Parameters: ! ! on entry ! ! a complex(lda, n) ! the hermitian matrix to be factored. only the ! diagonal and upper triangle are used. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! on return ! ! a 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. ! ! rcond real ! 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. if info /= 0 , rcond is unchanged. ! ! z complex(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) . ! if info /= 0 , z is unchanged. ! ! info integer ! = 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 jm1 integer k integer kb integer kp1 real rcond real s real scasum real sm complex t complex wk complex wkm 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),1),0.0e0) jm1 = j - 1 do i = 1, jm1 z(i) = cmplx(real(z(i))+cabs1(a(i,j)),0.0e0) end do end do anorm = 0.0e0 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 = (1.0e0,0.0e0) z(1:n) = (0.0e0,0.0e0) do k = 1, n if (cabs1(z(k)) /= 0.0e0) then ek = csign1(ek,-z(k)) end if if ( cabs1(ek-z(k)) > real(a(k,k)) ) then s = real(a(k,k))/cabs1(ek-z(k)) call csscal(n,s,z,1) ek = cmplx(s,0.0e0)*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 do j = kp1, n z(j) = z(j) + t*conjg(a(k,j)) end do end if end if z(k) = wk end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ! ! Solve r*y = w ! do kb = 1, n k = n + 1 - kb if (cabs1(z(k)) > real(a(k,k))) then s = real(a(k,k))/cabs1(z(k)) call csscal(n,s,z,1) end if z(k) = z(k)/a(k,k) t = -z(k) call caxpy(k-1,t,a(1,k),1,z(1),1) end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = 1.0e0 ! ! Solve hermitian(r)*v = y ! do k = 1, n z(k) = z(k) - cdotc(k-1,a(1,k),1,z(1),1) if ( cabs1(z(k)) > real(a(k,k)) ) then s = real(a(k,k))/cabs1(z(k)) call csscal(n,s,z,1) ynorm = s*ynorm end if z(k) = z(k)/a(k,k) end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = s*ynorm ! ! Solve r*z = v ! do kb = 1, n k = n + 1 - kb if ( cabs1(z(k)) > real(a(k,k)) ) then s = real(a(k,k))/cabs1(z(k)) call csscal(n,s,z,1) ynorm = s*ynorm end if z(k) = z(k)/a(k,k) t = -z(k) call caxpy(k-1,t,a(1,k),1,z(1),1) end do ! ! Make ZNORM = 1. ! s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = s*ynorm if ( anorm /= 0.0e0 ) then rcond = ynorm/anorm else rcond = 0.0e0 end if return end subroutine cpodi ( a, lda, n, det, job ) ! !******************************************************************************* ! !! CPODI computes the determinant and inverse of a certain ! complex hermitian positive definite matrix (see below) ! using the factors computed by cpoco, cpofa or cqrdc. ! ! error condition ! ! 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 . ! ! Parameters: ! ! on entry ! ! a complex(lda, n) ! the output a from cpoco or cpofa ! or the output x from cqrdc. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! job integer ! = 11 both determinant and inverse. ! = 01 inverse only. ! = 10 determinant only. ! ! on return ! ! a 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. ! ! det real(2) ! determinant of a or of hermitian(x)*x if requested. ! otherwise not referenced. ! 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 a(lda,n) real det(2) integer i integer j integer jm1 integer job integer k integer kp1 real s complex t ! ! Compute determinant ! if ( job/10 /= 0) then det(1) = 1.0e0 det(2) = 0.0e0 s = 10.0e0 do i = 1, n det(1) = real(a(i,i))**2*det(1) if (det(1) == 0.0e0) then exit end if do while (det(1) < 1.0e0) det(1) = s*det(1) det(2) = det(2) - 1.0e0 end do do while ( det(1) >= s ) det(1) = det(1)/s det(2) = det(2) + 1.0e0 end do end do end if ! ! Compute inverse(r) ! if ( mod(job,10) /= 0 ) then do k = 1, n a(k,k) = (1.0e0,0.0e0)/a(k,k) t = -a(k,k) call cscal(k-1,t,a(1,k),1) do j = k+1, n t = a(k,j) a(k,j) = (0.0e0,0.0e0) call caxpy(k,t,a(1,k),1,a(1,j),1) end do end do ! ! Form inverse(r) * hermitian(inverse(r)) ! do j = 1, n do k = 1, j-1 t = conjg(a(k,j)) call caxpy(k,t,a(1,j),1,a(1,k),1) end do t = conjg(a(j,j)) call cscal(j,t,a(1,j),1) end do end if return end subroutine cpofa ( a, lda, n, info ) ! !******************************************************************************* ! !! 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) . ! ! Parameters: ! ! on entry ! ! a complex(lda, n) ! the hermitian matrix to be factored. only the ! diagonal and upper triangle are used. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! on return ! ! a 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. ! ! info integer ! = 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.0e0 do k = 1, j-1 t = a(k,j) - cdotc(k-1,a(1,k),1,a(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.0e0 .or. aimag(a(j,j)) /= 0.0e0 ) then info = j exit end if a(j,j) = cmplx ( sqrt ( s ), 0.0e0 ) end do return end subroutine cposl ( a, lda, n, b ) ! !******************************************************************************* ! !! 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 /= 0) go to ... ! do j = 1, p ! call cposl(a,lda,n,c(1,j)) ! end do ! ! Parameters: ! ! on entry ! ! a complex(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 integer kb complex t ! ! Solve hermitian(r)*y = b ! do k = 1, n t = cdotc(k-1,a(1,k),1,b(1),1) b(k) = (b(k) - t)/a(k,k) end do ! ! Solve r*x = y ! do kb = 1, n k = n + 1 - kb b(k) = b(k)/a(k,k) t = -b(k) call caxpy(k-1,t,a(1,k),1,b(1),1) end do return end subroutine cppco ( ap, n, rcond, z, info ) ! !******************************************************************************* ! !! CPPCO factors a complex hermitian positive definite matrix ! stored in packed form ! and estimates the condition of the matrix. ! ! 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. ! ! 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 ! ! Parameters: ! ! on entry ! ! ap complex (n*(n+1)/2) ! 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 . ! see comments below for details. ! ! Input, integer N, the order of the matrix. ! ! on return ! ! ap an upper triangular matrix r , stored in packed ! form, so that a = hermitian(r)*r . ! if info /= 0 , the factorization is not complete. ! ! rcond real ! 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. if info /= 0 , rcond is unchanged. ! ! z complex(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. ! ! info integer ! = 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 jm1 integer k integer kb integer kj integer kk integer kp1 real rcond real s real scasum real sm complex t complex wk complex wkm real ynorm complex z(n) ! ! Find norm of A. ! j1 = 1 do j = 1, n z(j) = cmplx(scasum(j,ap(j1),1),0.0e0) ij = j1 j1 = j1 + j jm1 = j - 1 do i = 1, jm1 z(i) = cmplx(real(z(i))+cabs1(ap(ij)),0.0e0) 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 = (1.0e0,0.0e0) z(1:n) = (0.0e0,0.0e0) kk = 0 do k = 1, n kk = kk + k if (cabs1(z(k)) /= 0.0e0) ek = csign1(ek,-z(k)) if ( cabs1(ek-z(k)) > real(ap(kk)) ) then s = real(ap(kk))/cabs1(ek-z(k)) call csscal(n,s,z,1) ek = cmplx(s,0.0e0)*ek end if wk = ek - z(k) wkm = -ek - z(k) s = cabs1(wk) sm = cabs1(wkm) wk = wk/ap(kk) wkm = wkm/ap(kk) kp1 = k + 1 kj = kk + k if ( kp1 <= n ) then do j = kp1, n sm = sm + cabs1(z(j)+wkm*conjg(ap(kj))) z(j) = z(j) + wk*conjg(ap(kj)) s = s + cabs1(z(j)) kj = kj + j end do if ( s < sm ) then t = wkm - wk wk = wkm kj = kk + k do j = kp1, n z(j) = z(j) + t*conjg(ap(kj)) kj = kj + j end do end if end if z(k) = wk end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ! ! Solve r*y = W. ! do kb = 1, n k = n + 1 - kb if ( cabs1(z(k)) > real(ap(kk)) ) then s = real(ap(kk))/cabs1(z(k)) call csscal(n,s,z,1) end if z(k) = z(k)/ap(kk) kk = kk - k t = -z(k) call caxpy(k-1,t,ap(kk+1),1,z(1),1) end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = 1.0e0 ! ! Solve hermitian(r)*v = Y. ! do k = 1, n z(k) = z(k) - cdotc(k-1,ap(kk+1),1,z(1),1) kk = kk + k if ( cabs1(z(k)) > real(ap(kk)) ) then s = real(ap(kk))/cabs1(z(k)) call csscal(n,s,z,1) ynorm = s*ynorm end if z(k) = z(k)/ap(kk) end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = s*ynorm ! ! Solve r*z = V. ! do kb = 1, n k = n + 1 - kb if ( cabs1(z(k)) > real(ap(kk)) ) then s = real(ap(kk))/cabs1(z(k)) call csscal(n,s,z,1) ynorm = s*ynorm end if z(k) = z(k)/ap(kk) kk = kk - k t = -z(k) call caxpy(k-1,t,ap(kk+1),1,z(1),1) end do ! ! Make ZNORM = 1. ! s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = s*ynorm if ( anorm /= 0.0e0 ) then rcond = ynorm / anorm else rcond = 0.0e0 end if return end subroutine cppdi ( ap, n, det, job ) ! !******************************************************************************* ! !! CPPDI computes the determinant and inverse ! of a complex hermitian positive definite matrix ! using the factors computed by cppco or cppfa . ! ! 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 . ! ! Parameters: ! ! on entry ! ! ap complex (n*(n+1)/2) ! the output from cppco or cppfa. ! ! Input, integer N, the order of the matrix. ! ! job integer ! = 11 both determinant and inverse. ! = 01 inverse only. ! = 10 determinant only. ! ! on return ! ! ap the upper triangular half of the inverse . ! the strict lower triangle is unaltered. ! ! det real(2) ! determinant of original matrix if requested. ! otherwise not referenced. ! determinant = det(1) * 10.0**det(2) ! with 1.0 <= det(1) < 10.0 ! or det(1) == 0.0 . ! implicit none ! integer n ! complex ap((n*(n+1))/2) real det(2) integer i integer ii integer j integer j1 integer jj integer jm1 integer job integer k integer k1 integer kj integer kk integer kp1 real s complex t ! ! Compute determinant. ! if (job/10 /= 0) then det(1) = 1.0e0 det(2) = 0.0e0 s = 10.0e0 ii = 0 do i = 1, n ii = ii + i det(1) = real(ap(ii))**2*det(1) if ( det(1) == 0.0e0 ) then exit end if do while ( det(1) < 1.0e0 ) det(1) = s*det(1) det(2) = det(2) - 1.0e0 end do do while ( det(1) >= s ) det(1) = det(1)/s det(2) = det(2) + 1.0e0 end do end do end if ! ! Compute inverse ( R ). ! if ( mod(job,10) /= 0 ) then kk = 0 do k = 1, n k1 = kk + 1 kk = kk + k ap(kk) = (1.0e0,0.0e0)/ap(kk) t = -ap(kk) call cscal(k-1,t,ap(k1),1) kp1 = k + 1 j1 = kk + 1 kj = kk + k do j = kp1, n t = ap(kj) ap(kj) = (0.0e0,0.0e0) call caxpy(k,t,ap(k1),1,ap(j1),1) j1 = j1 + j kj = kj + j end do end do ! ! Form inverse ( R ) * hermitian ( inverse ( R ) ). ! jj = 0 do j = 1, n j1 = jj + 1 jj = jj + j jm1 = j - 1 k1 = 1 kj = j1 do k = 1, jm1 t = conjg(ap(kj)) call caxpy(k,t,ap(j1),1,ap(k1),1) k1 = k1 + k kj = kj + 1 end do t = conjg(ap(jj)) call cscal(j,t,ap(j1),1) end do end if return end subroutine cppfa ( ap, n, info ) ! !******************************************************************************* ! !! CPPFA factors a complex hermitian positive definite packed matrix. ! ! ! Discussion: ! ! 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 ! ! Parameters: ! ! on entry ! ! ap complex (n*(n+1)/2) ! 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 . ! see comments below for details. ! ! Input, integer N, the order of the matrix. ! ! on return ! ! ap an upper triangular matrix r , stored in packed ! form, so that a = hermitian(r)*r . ! ! info integer ! = 0 for normal return. ! = k if the leading minor of order k is not ! positive definite. ! implicit none ! integer n ! complex ap((n*(n+1))/2) complex cdotc integer info integer j integer jj integer jm1 integer k integer kj integer kk real s complex t ! info = 0 jj = 0 do j = 1, n s = 0.0e0 jm1 = j - 1 kj = jj kk = 0 do k = 1, jm1 kj = kj + 1 t = ap(kj) - cdotc(k-1,ap(kk+1),1,ap(jj+1),1) kk = kk + k t = t/ap(kk) ap(kj) = t s = s + real(t*conjg(t)) end do jj = jj + j s = real(ap(jj)) - s if ( s <= 0.0e0 .or. aimag(ap(jj) ) /= 0.0e0 ) then info = j exit end if ap(jj) = cmplx(sqrt(s),0.0e0) end do return end subroutine cppsl ( ap, n, b ) ! !******************************************************************************* ! !! CPPSL solves the complex hermitian positive definite system ! a * x = b ! using the factors computed by cppco or cppfa. ! ! error condition ! ! 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 cppco(ap,n,rcond,z,info) ! if (rcond is too small .or. info /= 0) go to ... ! do j = 1, p ! call cppsl(ap,n,c(1,j)) ! end do ! ! Parameters: ! ! on entry ! ! ap complex (n*(n+1)/2) ! the output from cppco or cppfa. ! ! 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 n ! complex ap((n*(n+1))/2) complex b(1) complex cdotc integer k integer kb integer kk complex t ! kk = 0 do k = 1, n t = cdotc(k-1,ap(kk+1),1,b(1),1) kk = kk + k b(k) = (b(k) - t)/ap(kk) end do do kb = 1, n k = n + 1 - kb b(k) = b(k)/ap(kk) kk = kk - k t = -b(k) call caxpy(k-1,t,ap(kk+1),1,b(1),1) end do return end subroutine cptsl ( n, d, e, b ) ! !******************************************************************************* ! !! CPTSL solves a positive definite tridiagonal linear system. ! ! ! Discussion; ! ! The system does not have to be factored first. ! ! Parameters: ! ! Input, integer N, the order of the matrix. ! ! Input/output, complex D(N). On input, the diagonal of the ! matrix. On output, this has been overwritten by other information. ! ! Input/output, complex E(N). On input, the superdiagonal entries of the ! matrix in locations E(1:N-1). On output, this has been overwritten ! by other information. ! ! Input/output, complex B(N). On input, the right hand side. ! On output, the solution. ! implicit none ! integer n ! complex b(n) complex d(n) complex e(n) integer k integer kbm1 integer ke integer kf integer kp1 integer nm1 integer nm1d2 complex t1 complex t2 ! ! Check for 1 x 1 case ! if ( n == 1 ) then b(1) = b(1) / d(1) return end if nm1 = n - 1 nm1d2 = ( n - 1 ) / 2 if ( n /= 2 ) then kbm1 = n - 1 ! ! Zero top half of subdiagonal and bottom half of superdiagonal. ! do k = 1, nm1d2 t1 = conjg(e(k))/d(k) d(k+1) = d(k+1) - t1*e(k) b(k+1) = b(k+1) - t1*b(k) t2 = e(kbm1)/d(kbm1+1) d(kbm1) = d(kbm1) - t2*conjg(e(kbm1)) b(kbm1) = b(kbm1) - t2*b(kbm1+1) kbm1 = kbm1 - 1 end do end if kp1 = nm1d2 + 1 ! ! Clean up for possible 2 x 2 block at center. ! if ( mod(n,2) == 0 ) then t1 = conjg(e(kp1))/d(kp1) d(kp1+1) = d(kp1+1) - t1*e(kp1) b(kp1+1) = b(kp1+1) - t1*b(kp1) kp1 = kp1 + 1 end if ! ! Back solve starting at the center, going towards the top and bottom. ! b(kp1) = b(kp1)/d(kp1) if ( n /= 2 ) then k = kp1 - 1 ke = kp1 + nm1d2 - 1 do kf = kp1, ke b(k) = (b(k) - e(k)*b(k+1))/d(k) b(kf+1) = (b(kf+1) - conjg(e(kf))*b(kf))/d(kf+1) k = k - 1 end do end if if ( mod ( n, 2 ) == 0 ) then b(1) = ( b(1) - e(1) * b(2) ) / d(1) end if return end subroutine cqrdc ( x, ldx, n, p, qraux, ipvt, work, job ) ! !******************************************************************************* ! !! CQRDC computes the QR factorization of an N by P complex matrix. ! ! ! Discussion: ! ! CQRDC uses Householder transformations to compute the QR factorization ! of an N by P matrix X. Column pivoting based on the 2-norms of the ! reduced columns may be performed at the user's option. ! ! Parameters: ! ! on entry ! ! x complex(ldx,p), where ldx >= n. ! x contains the matrix whose decomposition is to be ! computed. ! ! Input, integer LDX, the leading dimension of X. ! ! Input, integer N, the number of rows of the matrix. ! ! p integer. ! p is the number of columns of the matrix x. ! ! ipvt integer(p). ! ipvt contains integers that control the selection ! of the pivot columns. the k-th column x(k) of x ! is placed in one of three classes according to the ! value of ipvt(k). ! ! if ipvt(k) > 0, then x(k) is an initial ! column. ! ! if ipvt(k) == 0, then x(k) is a free column. ! ! if ipvt(k) < 0, then x(k) is a final column. ! ! before the decomposition is computed, initial columns ! are moved to the beginning of the array x and final ! columns to the end. both initial and final columns ! are frozen in place during the computation and only ! free columns are moved. at the k-th stage of the ! reduction, if x(k) is occupied by a free column ! it is interchanged with the free column of largest ! reduced norm. ipvt is not referenced if ! job == 0. ! ! work complex(p). ! work is a work array. work is not referenced if ! job == 0. ! ! job integer. ! job is an integer that initiates column pivoting. ! if job == 0, no pivoting is done. ! if job /= 0, pivoting is done. ! ! on return ! ! x x contains in its upper triangle the upper ! triangular matrix r of the qr factorization. ! below its diagonal x contains information from ! which the unitary part of the decomposition ! can be recovered. note that if pivoting has ! been requested, the decomposition is not that ! of the original matrix x but that of x ! with its columns permuted as described by ipvt. ! ! qraux complex(p). ! qraux contains further information required to recover ! the unitary part of the decomposition. ! ! ipvt ipvt(k) contains the index of the column of the ! original matrix that has been interchanged into ! the k-th column, if pivoting was requested. ! implicit none ! integer ldx integer n integer p ! real cabs1 complex cdotc complex csign integer ipvt(p) integer j integer jj integer job integer jp integer l integer lp1 integer lup integer maxj real maxnrm logical negj complex nrmxl integer pl integer pu complex qraux(p) real scnrm2 logical swapj complex t real tt complex work(p) complex x(ldx,p) ! pl = 1 pu = 0 if ( job /= 0 ) then ! ! Pivoting has been requested. Rearrange the columns according to IPVT. ! do j = 1, p swapj = ipvt(j) > 0 negj = ipvt(j) < 0 ipvt(j) = j if (negj) ipvt(j) = -j if ( swapj ) then if (j /= pl) call cswap(n,x(1,pl),1,x(1,j),1) ipvt(j) = ipvt(pl) ipvt(pl) = j pl = pl + 1 end if end do pu = p do jj = 1, p j = p - jj + 1 if ( ipvt(j) < 0 ) then ipvt(j) = -ipvt(j) if ( j /= pu ) then call cswap(n,x(1,pu),1,x(1,j),1) jp = ipvt(pu) ipvt(pu) = ipvt(j) ipvt(j) = jp end if pu = pu - 1 end if end do end if ! ! Compute the norms of the free columns. ! do j = pl, pu qraux(j) = cmplx(scnrm2(n,x(1,j),1),0.0e0) work(j) = qraux(j) end do ! ! Perform the Householder reduction of X. ! lup = min(n,p) do l = 1, lup ! ! locate the column of largest norm and bring it ! into the pivot position. ! if ( l >= pl .and. l < pu ) then maxnrm = 0.0e0 maxj = l do j = l, pu if ( real(qraux(j)) > maxnrm) then maxnrm = real(qraux(j)) maxj = j end if end do if ( maxj /= l ) then call cswap(n,x(1,l),1,x(1,maxj),1) qraux(maxj) = qraux(l) work(maxj) = work(l) jp = ipvt(maxj) ipvt(maxj) = ipvt(l) ipvt(l) = jp end if end if qraux(l) = (0.0e0,0.0e0) if ( l /= n ) then ! ! Compute the Householder transformation for column L. ! nrmxl = cmplx(scnrm2(n-l+1,x(l,l),1),0.0e0) if ( cabs1(nrmxl) /= 0.0e0 ) then if (cabs1(x(l,l)) /= 0.0e0) then nrmxl = csign(nrmxl,x(l,l)) end if call cscal(n-l+1,(1.0e0,0.0e0)/nrmxl,x(l,l),1) x(l,l) = (1.0e0,0.0e0) + x(l,l) ! ! Apply the transformation to the remaining columns, ! updating the norms. ! lp1 = l + 1 do j = lp1, p t = -cdotc(n-l+1,x(l,l),1,x(l,j),1)/x(l,l) call caxpy(n-l+1,t,x(l,l),1,x(l,j),1) if ( j < pl .or. j > pu ) then cycle end if if ( cabs1(qraux(j)) == 0.0e0 ) then cycle end if tt = 1.0e0 - (cabs(x(l,j))/real(qraux(j)))**2 tt = max ( tt, 0.0e0 ) t = cmplx(tt,0.0e0) tt = 1.0e0 + 0.05e0*tt*(real(qraux(j))/real(work(j)))**2 if ( tt /= 1.0e0 ) then qraux(j) = qraux(j)*csqrt(t) else qraux(j) = cmplx(scnrm2(n-l,x(l+1,j),1),0.0e0) work(j) = qraux(j) end if end do ! ! Save the transformation. ! qraux(l) = x(l,l) x(l,l) = -nrmxl end if end if end do return end subroutine cqrsl ( x, ldx, n, k, qraux, y, qy, qty, b, rsd, xb, job, info ) ! !******************************************************************************* ! !! CQRSL applies the output of cqrdc to compute coordinate ! transformations, projections, and least squares solutions. ! for k <= min(n,p), let xk be the matrix ! ! xk = (x(ipvt(1)),x(ipvt(2)), ... ,x(ipvt(k))) ! ! formed from columnns ipvt(1), ... ,ipvt(k) of the original ! n x p matrix x that was input to cqrdc (if no pivoting was ! done, xk consists of the first k columns of x in their ! original order). cqrdc produces a factored unitary matrix q ! and an upper triangular matrix r such that ! ! xk = q * (r) ! (0) ! ! this information is contained in coded form in the arrays ! x and qraux. ! ! Parameters: ! ! on entry ! ! x complex(ldx,p). ! x contains the output of cqrdc. ! ! Input, integer LDX, the leading dimension of X. ! ! Input, integer N, the number of rows of the matrix XK, which ! must have the same value as N in CQRDC. ! ! k integer. ! k is the number of columns of the matrix xk. k ! must nnot be greater than min(n,p), where p is the ! same as in the calling sequence to cqrdc. ! ! qraux complex(p). ! qraux contains the auxiliary output from cqrdc. ! ! y complex(n) ! y contains an n-vector that is to be manipulated ! by cqrsl. ! ! job integer. ! job specifies what is to be computed. job has ! the decimal expansion abcde, with the following ! meaning. ! ! if a/=0, compute qy. ! if b,c,d, or e /= 0, compute qty. ! if c/=0, compute b. ! if d/=0, compute rsd. ! if e/=0, compute xb. ! ! note that a request to compute b, rsd, or xb ! automatically triggers the computation of qty, for ! which an array must be provided in the calling ! sequence. ! ! on return ! ! qy complex(n). ! qy conntains q*y, if its computation has been ! requested. ! ! qty complex(n). ! qty contains hermitian(q)*y, if its computation has ! been requested. here hermitian(q) is the conjugate ! transpose of the matrix q. ! ! b complex(k) ! b contains the solution of the least squares problem ! ! minimize norm2(y - xk*b), ! ! if its computation has been requested. (note that ! if pivoting was requested in cqrdc, the j-th ! component of b will be associated with column ipvt(j) ! of the original matrix x that was input into cqrdc.) ! ! rsd complex(n). ! rsd contains the least squares residual y - xk*b, ! if its computation has been requested. rsd is ! also the orthogonal projection of y onto the ! orthogonal complement of the column space of xk. ! ! xb complex(n). ! xb contains the least squares approximation xk*b, ! if its computation has been requested. xb is also ! the orthogonal projection of y onto the column space ! of x. ! ! info integer. ! info is zero unless the computation of b has ! been requested and r is exactly singular. in ! this case, info is the index of the first zero ! diagonal element of r and b is left unaltered. ! ! the parameters qy, qty, b, rsd, and xb are not referenced ! if their computation is not requested and in this case ! can be replaced by dummy variables in the calling program. ! to save storage, the user may in some cases use the same ! array for different parameters in the calling sequence. a ! frequently occuring example is when one wishes to compute ! any of b, rsd, or xb and does not need y or qty. in this ! case one may identify y, qty, and one of b, rsd, or xb, while ! providing separate arrays for anything else that is to be ! computed. thus the calling sequence ! ! call cqrsl(x,ldx,n,k,qraux,y,dum,y,b,y,dum,110,info) ! ! will result in the computation of b and rsd, with rsd ! overwriting y. more generally, each item in the following ! list contains groups of permissible identifications for ! a single callinng sequence. ! ! 1. (y,qty,b) (rsd) (xb) (qy) ! ! 2. (y,qty,rsd) (b) (xb) (qy) ! ! 3. (y,qty,xb) (b) (rsd) (qy) ! ! 4. (y,qy) (qty,b) (rsd) (xb) ! ! 5. (y,qy) (qty,rsd) (b) (xb) ! ! 6. (y,qy) (qty,xb) (b) (rsd) ! ! in any group the value returned in the array allocated to ! the group corresponds to the last member of the group. ! implicit none ! integer k integer ldx integer n ! complex b(k) real cabs1 logical cb complex cdotc logical cqty logical cqy logical cr logical cxb integer i integer info integer j integer jj integer job integer ju integer kp1 complex qraux(*) complex qty(n) complex qy(n) complex rsd(n) complex t complex temp complex x(ldx,*) complex xb(n) complex y(n) ! info = 0 ! ! Determine what is to be computed. ! cqy = job/10000 /= 0 cqty = mod(job,10000) /= 0 cb = mod(job,1000)/100 /= 0 cr = mod(job,100)/10 /= 0 cxb = mod(job,10) /= 0 ju = min(k,n-1) ! ! Special action when N=1. ! if ( ju == 0 ) then if (cqy) qy(1) = y(1) if (cqty) qty(1) = y(1) if (cxb) xb(1) = y(1) if ( cb ) then if ( cabs1(x(1,1)) == 0.0e0 ) then info = 1 else b(1) = y(1)/x(1,1) end if end if if (cr) rsd(1) = (0.0e0,0.0e0) return end if ! ! Set up to compute QY or QTY. ! if (cqy) call ccopy(n,y,1,qy,1) if (cqty) call ccopy(n,y,1,qty,1) ! ! Compute QY. ! if (cqy) then do jj = 1, ju j = ju - jj + 1 if ( cabs1(qraux(j)) /= 0.0e0 ) then temp = x(j,j) x(j,j) = qraux(j) t = -cdotc(n-j+1,x(j,j),1,qy(j),1)/x(j,j) call caxpy(n-j+1,t,x(j,j),1,qy(j),1) x(j,j) = temp end if end do end if ! ! Compute hermitian ( A ) * Y. ! if ( cqty ) then do j = 1, ju if ( cabs1(qraux(j)) /= 0.0e0 ) then temp = x(j,j) x(j,j) = qraux(j) t = -cdotc(n-j+1,x(j,j),1,qty(j),1)/x(j,j) call caxpy(n-j+1,t,x(j,j),1,qty(j),1) x(j,j) = temp end if end do end if ! ! Set up to compute B, RSD, or XB. ! if (cb) call ccopy(k,qty,1,b,1) kp1 = k + 1 if (cxb) call ccopy(k,qty,1,xb,1) if (cr .and. k < n) call ccopy(n-k,qty(kp1),1,rsd(kp1),1) if ( cxb ) then xb(kp1:n) = (0.0e0,0.0e0) end if if (cr) then rsd(1:k) = (0.0e0,0.0e0) end if ! ! Compute B. ! if ( cb ) then do jj = 1, k j = k - jj + 1 if ( cabs1(x(j,j)) == 0.0e0 ) then info = j exit end if b(j) = b(j) / x(j,j) if ( j /= 1 ) then t = -b(j) call caxpy(j-1,t,x(1,j),1,b,1) end if end do end if if ( cr .or. cxb ) then ! ! Compute RSD or XB as required. ! do jj = 1, ju j = ju - jj + 1 if ( cabs1(qraux(j)) /= 0.0e0 ) then temp = x(j,j) x(j,j) = qraux(j) if ( cr ) then t = -cdotc(n-j+1,x(j,j),1,rsd(j),1)/x(j,j) call caxpy(n-j+1,t,x(j,j),1,rsd(j),1) end if if ( cxb ) then t = -cdotc(n-j+1,x(j,j),1,xb(j),1)/x(j,j) call caxpy(n-j+1,t,x(j,j),1,xb(j),1) end if x(j,j) = temp end if end do end if return end subroutine csico ( a, lda, n, ipvt, rcond, z ) ! !******************************************************************************* ! !! CSICO factors a complex symmetric matrix by elimination with ! symmetric pivoting and estimates the condition of the matrix. ! ! if rcond is not needed, csifa is slightly faster. ! to solve a*x = b , follow csico by csisl. ! to compute inverse(a)*c , follow csico by csisl. ! to compute inverse(a) , follow csico by csidi. ! to compute determinant(a) , follow csico by csidi. ! ! Parameters: ! ! on entry ! ! a complex(lda, n) ! the symmetric matrix to be factored. ! only the diagonal and upper triangle are used. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! on return ! ! a a block diagonal matrix and the multipliers which ! were used to obtain it. ! the factorization can be written a = u*d*trans(u) ! where u is a product of permutation and unit ! upper triangular matrices , trans(u) is the ! transpose of u , and d is block diagonal ! with 1 by 1 and 2 by 2 blocks. ! ! ipvt integer(n) ! an integer vector of pivot indices. ! ! rcond real ! 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. ! ! z complex(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 cdotu complex csign1 complex denom complex ek integer i integer info integer ipvt(n) integer j integer jm1 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),1),0.0e0) jm1 = j - 1 do i = 1, jm1 z(i) = cmplx(real(z(i))+cabs1(a(i,j)),0.0e0) end do end do anorm = 0.0e0 do j = 1, n anorm = max (anorm,real(z(j))) end do ! ! Factor. ! call csifa(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 = (1.0e0,0.0e0) z(1:n) = (0.0e0,0.0e0) k = n do while ( k > 0 ) ks = 1 if (ipvt(k) < 0) ks = 2 kp = abs(ipvt(k)) kps = k + 1 - ks if ( kp /= kps ) then call c_swap ( z(kps), z(kp) ) end if if ( cabs1(z(k)) /= 0.0e0 ) then ek = csign1(ek,z(k)) end if z(k) = z(k) + ek call caxpy(k-ks,z(k),a(1,k),1,z(1),1) if ( ks /= 1 ) then if (cabs1(z(k-1)) /= 0.0e0) ek = csign1(ek,z(k-1)) z(k-1) = z(k-1) + ek call caxpy(k-ks,z(k-1),a(1,k-1),1,z(1),1) end if if ( ks /= 2 ) then if ( cabs1(z(k)) > cabs1(a(k,k)) ) then s = cabs1(a(k,k))/cabs1(z(k)) call csscal(n,s,z,1) ek = cmplx(s,0.0e0)*ek end if if (cabs1(a(k,k)) /= 0.0e0) z(k) = z(k)/a(k,k) if (cabs1(a(k,k)) == 0.0e0) z(k) = (1.0e0,0.0e0) else ak = a(k,k)/a(k-1,k) akm1 = a(k-1,k-1)/a(k-1,k) bk = z(k)/a(k-1,k) bkm1 = z(k-1)/a(k-1,k) denom = ak*akm1 - 1.0e0 z(k) = (akm1*bk - bkm1)/denom z(k-1) = (ak*bkm1 - bk)/denom end if k = k - ks end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ! ! Solve U' * Y = W. ! k = 1 do while ( k <= n ) ks = 1 if (ipvt(k) < 0) ks = 2 if ( k /= 1 ) then z(k) = z(k) + cdotu(k-1,a(1,k),1,z(1),1) if (ks == 2) then z(k+1) = z(k+1) + cdotu(k-1,a(1,k+1),1,z(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 k = k + ks end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = 1.0e0 ! ! Solve u*d*v = Y. ! k = n do while ( k > 0 ) ks = 1 if (ipvt(k) < 0) ks = 2 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 call caxpy(k-ks,z(k),a(1,k),1,z(1),1) if (ks == 2) call caxpy(k-ks,z(k-1),a(1,k-1),1,z(1),1) end if if ( ks /= 2 ) then if ( cabs1(z(k)) > cabs1(a(k,k)) ) then s = cabs1(a(k,k))/cabs1(z(k)) call csscal(n,s,z,1) ynorm = s*ynorm end if if (cabs1(a(k,k)) /= 0.0e0) z(k) = z(k)/a(k,k) if (cabs1(a(k,k)) == 0.0e0) z(k) = (1.0e0,0.0e0) else ak = a(k,k)/a(k-1,k) akm1 = a(k-1,k-1)/a(k-1,k) bk = z(k)/a(k-1,k) bkm1 = z(k-1)/a(k-1,k) denom = ak*akm1 - 1.0e0 z(k) = (akm1*bk - bkm1)/denom z(k-1) = (ak*bkm1 - bk)/denom end if k = k - ks end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = s*ynorm ! ! Solve U'*z = V. ! k = 1 do while ( k <= n ) ks = 1 if (ipvt(k) < 0) ks = 2 if (k /= 1) then z(k) = z(k) + cdotu(k-1,a(1,k),1,z(1),1) if (ks == 2) then z(k+1) = z(k+1) + cdotu(k-1,a(1,k+1),1,z(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 k = k + ks end do ! ! Make ZNORM = 1. ! s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = s*ynorm if ( anorm /= 0.0e0 ) then rcond = ynorm / anorm else rcond = 0.0e0 end if return end subroutine csidi ( a, lda, n, ipvt, det, work, job ) ! !******************************************************************************* ! !! CSIDI computes the determinant and inverse ! of a complex symmetric matrix using the factors from csifa. ! ! Discussion: ! ! a division by zero may occur if the inverse is requested ! and csico has set rcond == 0.0 ! or csifa has set info /= 0 . ! ! Parameters: ! ! on entry ! ! a complex(lda,n) ! the output from csifa. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! ipvt integer(n) ! the pivot vector from csifa. ! ! Workspace, complex WORK(N). ! ! job integer ! job has the decimal expansion ab where ! if b /= 0, the inverse is computed, ! if a /= 0, the determinant is computed, ! ! for example, job = 11 gives both. ! ! on return ! ! variables not requested by job are not used. ! ! a contains the upper triangle of the inverse of ! the original matrix. the strict lower triangle ! is never referenced. ! ! det complex(2) ! determinant of original matrix. ! determinant = det(1) * 10.0**det(2) ! with 1.0 <= abs(det(1)) < 10.0 ! or det(1) = 0.0. ! implicit none ! integer lda integer n ! complex a(lda,n) complex ak complex akkp1 complex akp1 real cabs1 complex cdotu complex d complex det(2) integer ipvt(n) integer j integer jb integer job integer k integer km1 integer ks integer kstep logical nodet logical noinv complex t complex temp real, parameter :: ten = 10.0E+00 complex work(n) ! noinv = mod(job,10) == 0 nodet = mod(job,100)/10 == 0 if ( .not. nodet ) then det(1) = (1.0e0,0.0e0) det(2) = (0.0e0,0.0e0) t = (0.0e0,0.0e0) do k = 1, n d = a(k,k) ! ! 2 by 2 block ! use det (d t) = (d/t * c - t) * t ! (t c) ! to avoid underflow/overflow troubles. ! take two passes through scaling. use t for flag. ! if (ipvt(k) <= 0 ) then if ( cabs1(t) == 0.0e0) then t = a(k,k+1) d = (d/t)*a(k+1,k+1) - t else d = t t = (0.0e0,0.0e0) end if end if det(1) = d*det(1) if ( cabs1(det(1)) /= 0.0e0 ) then do while ( cabs1(det(1)) < 1.0e0 ) det(1) = cmplx(ten,0.0e0)*det(1) det(2) = det(2) - (1.0e0,0.0e0) end do do while ( cabs1(det(1)) >= ten ) det(1) = det(1)/cmplx(ten,0.0e0) det(2) = det(2) + (1.0e0,0.0e0) end do end if end do end if ! ! Compute inverse ( A ). ! if ( .not. noinv ) then k = 1 do while ( k <= n ) km1 = k - 1 ! ! 1 by 1 ! if ( ipvt(k) >= 0) then a(k,k) = (1.0e0,0.0e0)/a(k,k) if ( km1 >= 1 ) then call ccopy(km1,a(1,k),1,work,1) do j = 1, km1 a(j,k) = cdotu(j,a(1,j),1,work,1) call caxpy(j-1,work(j),a(1,j),1,a(1,k),1) end do a(k,k) = a(k,k) + cdotu(km1,work,1,a(1,k),1) end if kstep = 1 ! ! 2 by 2 ! else t = a(k,k+1) ak = a(k,k)/t akp1 = a(k+1,k+1)/t akkp1 = a(k,k+1)/t d = t*(ak*akp1 - (1.0e0,0.0e0)) a(k,k) = akp1/d a(k+1,k+1) = ak/d a(k,k+1) = -akkp1/d if (km1 >= 1) then call ccopy(km1,a(1,k+1),1,work,1) do j = 1, km1 a(j,k+1) = cdotu(j,a(1,j),1,work,1) call caxpy(j-1,work(j),a(1,j),1,a(1,k+1),1) end do a(k+1,k+1) = a(k+1,k+1) + cdotu(km1,work,1,a(1,k+1),1) a(k,k+1) = a(k,k+1) + cdotu(km1,a(1,k),1,a(1,k+1),1) call ccopy(km1,a(1,k),1,work,1) do j = 1, km1 a(j,k) = cdotu(j,a(1,j),1,work,1) call caxpy(j-1,work(j),a(1,j),1,a(1,k),1) end do a(k,k) = a(k,k) + cdotu(km1,work,1,a(1,k),1) end if kstep = 2 end if ! ! Swap. ! ks = abs(ipvt(k)) if ( ks /= k ) then call cswap(ks,a(1,ks),1,a(1,k),1) do jb = ks, k j = k + ks - jb temp = a(j,k) a(j,k) = 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 csifa ( a, lda, n, ipvt, info ) ! !******************************************************************************* ! !! CSIFA factors a complex symmetric matrix ! by elimination with symmetric pivoting. ! ! to solve a*x = b , follow csifa by csisl. ! to compute inverse(a)*c , follow csifa by csisl. ! to compute determinant(a) , follow csifa by csidi. ! to compute inverse(a) , follow csifa by csidi. ! ! Parameters: ! ! on entry ! ! a complex(lda,n) ! the symmetric matrix to be factored. ! only the diagonal and upper triangle are used. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! on return ! ! a a block diagonal matrix and the multipliers which ! were used to obtain it. ! the factorization can be written a = u*d*trans(u) ! where u is a product of permutation and unit ! upper triangular matrices , trans(u) is the ! transpose of u , and d is block diagonal ! with 1 by 1 and 2 by 2 blocks. ! ! ipvt integer(n) ! an integer vector of pivot indices. ! ! info integer ! = 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 csisl or csidi 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 imaxp1 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.0e0 + sqrt(17.0e0))/8.0e0 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.0e0) info = 1 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) colmax = cabs1(a(imax,k)) if ( absakk >= alpha*colmax) then kstep = 1 swap = .false. ! ! Determine the largest off-diagonal element in row IMAX. ! else rowmax = 0.0e0 imaxp1 = imax + 1 do j = imaxp1, k rowmax = max (rowmax,cabs1(a(imax,j))) end do if ( imax /= 1 ) then jmax = icamax(imax-1,a(1,imax),1) rowmax = max (rowmax,cabs1(a(jmax,imax))) end if if ( cabs1(a(imax,imax)) >= alpha*rowmax ) then kstep = 1 swap = .true. else if ( absakk >= alpha*colmax*(colmax/rowmax) ) 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.0e0) then ipvt(k) = k info = k k = k - kstep cycle end if if ( kstep /= 2 ) then ! ! 1 x 1 pivot block. ! ! Perform an interchange. ! if ( swap ) then call cswap(imax,a(1,imax),1,a(1,k),1) do jj = imax, k j = k + imax - jj t = a(j,k) a(j,k) = 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 = mulk call caxpy(j,t,a(1,k),1,a(1,j),1) a(j,k) = mulk end do ! ! Set the pivot array. ! ipvt(k) = k if (swap) ipvt(k) = imax ! ! 2 x 2 pivot block. ! else if ( swap ) then call cswap(imax,a(1,imax),1,a(1,k-1),1) do jj = imax, km1 j = km1 + imax - jj t = a(j,k-1) a(j,k-1) = 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 (km2 /= 0) then ak = a(k,k)/a(k-1,k) akm1 = a(k-1,k-1)/a(k-1,k) denom = 1.0e0 - ak*akm1 do jj = 1, km2 j = km1 - jj bk = a(j,k)/a(k-1,k) bkm1 = a(j,k-1)/a(k-1,k) mulk = (akm1*bk - bkm1)/denom mulkm1 = (ak*bkm1 - bk)/denom t = mulk call caxpy(j,t,a(1,k),1,a(1,j),1) t = mulkm1 call caxpy(j,t,a(1,k-1),1,a(1,j),1) a(j,k) = mulk a(j,k-1) = mulkm1 end do end if ! ! Set the pivot array. ! ipvt(k) = 1 - k if (swap) ipvt(k) = -imax ipvt(k-1) = ipvt(k) end if k = k - kstep end do return end subroutine csisl ( a, lda, n, ipvt, b ) ! !******************************************************************************* ! !! CSISL solves a complex symmetric system that was factored by CSIFA. ! ! ! Discussion: ! ! A division by zero may occur if csico has set rcond == 0.0 ! or csifa has set info /= 0 . ! ! To compute inverse(a) * c where c is a matrix with p columns ! ! call csifa(a,lda,n,ipvt,info) ! if (info /= 0) go to ... ! do j = 1, p ! call csisl(a,lda,n,ipvt,c(1,j)) ! end do ! ! Parameters: ! ! Input, complex A(LDA,N), the output from CSICO or CSIFA. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! Input, integer IPVT(N), the pivot vector from CSICO or CSIFA. ! ! 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 cdotu complex denom integer ipvt(n) integer k integer kp complex temp ! ! Loop backward applying the transformations and d inverse to b. ! k = n do while ( k > 0 ) ! ! 1 x 1 pivot block. ! if (ipvt(k) >= 0) then if ( k /= 1 ) then kp = ipvt(k) if ( kp /= k ) then temp = b(k) b(k) = b(kp) b(kp) = temp end if call caxpy(k-1,b(k),a(1,k),1,b(1),1) end if 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 temp = b(k-1) b(k-1) = b(kp) b(kp) = temp end if call caxpy(k-2,b(k),a(1,k),1,b(1),1) call caxpy(k-2,b(k-1),a(1,k-1),1,b(1),1) end if ak = a(k,k)/a(k-1,k) akm1 = a(k-1,k-1)/a(k-1,k) bk = b(k)/a(k-1,k) bkm1 = b(k-1)/a(k-1,k) denom = ak*akm1 - 1.0e0 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 ) if ( ipvt(k) >= 0 ) then ! ! 1 x 1 pivot block. ! if ( k /= 1 ) then b(k) = b(k) + cdotu(k-1,a(1,k),1,b(1),1) kp = ipvt(k) if ( kp /= k ) then call c_swap ( b(k), b(kp) ) end if end if k = k + 1 ! ! 2 x 2 pivot block. ! else if (k /= 1) then b(k) = b(k) + cdotu(k-1,a(1,k),1,b(1),1) b(k+1) = b(k+1) + cdotu(k-1,a(1,k+1),1,b(1),1) kp = abs(ipvt(k)) if ( kp /= k ) then temp = b(k) b(k) = b(kp) b(kp) = temp end if end if k = k + 2 end if end do return end subroutine cspco ( ap, n, ipvt, rcond, z ) ! !******************************************************************************* ! !! CSPCO factors a complex symmetric matrix stored in packed ! form by elimination with symmetric pivoting and estimates ! the condition of the matrix. ! ! if rcond is not needed, cspfa is slightly faster. ! to solve a*x = b , follow cspco by cspsl. ! to compute inverse(a)*c , follow cspco by cspsl. ! to compute inverse(a) , follow cspco by cspdi. ! to compute determinant(a) , follow cspco by cspdi. ! ! packed storage ! ! the following program segment will pack the upper ! triangle of a symmetric matrix. ! ! k = 0 ! do j = 1, n ! do i = 1, j ! k = k + 1 ! ap(k) = a(i,j) ! end do ! end do ! ! Parameters: ! ! on entry ! ! ap complex (n*(n+1)/2) ! the packed form of a symmetric matrix a . the ! columns of the upper triangle are stored sequentially ! in a one-dimensional array of length n*(n+1)/2 . ! see comments below for details. ! ! Input, integer N, the order of the matrix. ! ! on return ! ! ap 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*trans(u) ! where u is a product of permutation and unit ! upper triangular matrices , trans(u) is the ! transpose of u , and d is block diagonal ! with 1 by 1 and 2 by 2 blocks. ! ! ipvt integer(n) ! an integer vector of pivot indices. ! ! rcond real ! 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. ! ! z complex(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 cdotu 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 jm1 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),1),0.0e0) ij = j1 j1 = j1 + j jm1 = j - 1 do i = 1, jm1 z(i) = cmplx(real(z(i))+cabs1(ap(ij)),0.0e0) ij = ij + 1 end do end do anorm = 0.0e0 do j = 1, n anorm = max (anorm,real(z(j))) end do ! ! Factor. ! call cspfa(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 = (1.0e0,0.0e0) z(1:n) = (0.0e0,0.0e0) k = n ik = (n*(n - 1))/2 do while ( k > 0 ) kk = ik + k ikm1 = ik - (k - 1) ks = 1 if (ipvt(k) < 0) ks = 2 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.0e0) ek = csign1(ek,z(k)) z(k) = z(k) + ek call caxpy(k-ks,z(k),ap(ik+1),1,z(1),1) if ( ks /= 1 ) then if (cabs1(z(k-1)) /= 0.0e0) ek = csign1(ek,z(k-1)) z(k-1) = z(k-1) + ek call caxpy(k-ks,z(k-1),ap(ikm1+1),1,z(1),1) end if if ( ks /= 2 ) then if (cabs1(z(k)) > cabs1(ap(kk))) then s = cabs1(ap(kk))/cabs1(z(k)) call csscal(n,s,z,1) ek = cmplx(s,0.0e0)*ek end if if (cabs1(ap(kk)) /= 0.0e0) z(k) = z(k)/ap(kk) if (cabs1(ap(kk)) == 0.0e0) z(k) = (1.0e0,0.0e0) else km1k = ik + k - 1 km1km1 = ikm1 + k - 1 ak = ap(kk)/ap(km1k) akm1 = ap(km1km1)/ap(km1k) bk = z(k)/ap(km1k) bkm1 = z(k-1)/ap(km1k) denom = ak*akm1 - 1.0e0 z(k) = (akm1*bk - bkm1)/denom z(k-1) = (ak*bkm1 - bk)/denom end if k = k - ks ik = ik - k if (ks == 2) ik = ik - (k + 1) end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ! ! Solve trans(u)*y = W. ! k = 1 ik = 0 do while ( k <= n ) ks = 1 if (ipvt(k) < 0) ks = 2 if (k /= 1) then z(k) = z(k) + cdotu(k-1,ap(ik+1),1,z(1),1) ikp1 = ik + k if (ks == 2) then z(k+1) = z(k+1) + cdotu(k-1,ap(ikp1+1),1,z(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) ik = ik + (k + 1) k = k + ks end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = 1.0e0 ! ! Solve u*d*v = Y. ! k = n ik = n*(n - 1)/2 do while ( k > 0 ) kk = ik + k ikm1 = ik - (k - 1) ks = 1 if (ipvt(k) < 0) ks = 2 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 call caxpy(k-ks,z(k),ap(ik+1),1,z(1),1) if (ks == 2) call caxpy(k-ks,z(k-1),ap(ikm1+1),1,z(1),1) end if if ( ks /= 2 ) then if (cabs1(z(k)) > cabs1(ap(kk))) then s = cabs1(ap(kk))/cabs1(z(k)) call csscal(n,s,z,1) ynorm = s*ynorm end if if (cabs1(ap(kk)) /= 0.0e0) z(k) = z(k)/ap(kk) if (cabs1(ap(kk)) == 0.0e0) z(k) = (1.0e0,0.0e0) else km1k = ik + k - 1 km1km1 = ikm1 + k - 1 ak = ap(kk)/ap(km1k) akm1 = ap(km1km1)/ap(km1k) bk = z(k)/ap(km1k) bkm1 = z(k-1)/ap(km1k) denom = ak*akm1 - 1.0e0 z(k) = (akm1*bk - bkm1)/denom z(k-1) = (ak*bkm1 - bk)/denom end if k = k - ks ik = ik - k if (ks == 2) ik = ik - (k + 1) end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = s*ynorm ! ! Solve U' * Z = V. ! k = 1 ik = 0 do while ( k <= n ) ks = 1 if (ipvt(k) < 0) ks = 2 if (k /= 1) then z(k) = z(k) + cdotu(k-1,ap(ik+1),1,z(1),1) ikp1 = ik + k if (ks == 2) then z(k+1) = z(k+1) + cdotu(k-1,ap(ikp1+1),1,z(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) ik = ik + (k + 1) k = k + ks end do ! ! Make ZNORM = 1. ! s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = s*ynorm if ( anorm /= 0.0e0 ) then rcond = ynorm / anorm else rcond = 0.0e0 end if return end subroutine cspdi ( ap, n, ipvt, det, work, job ) ! !******************************************************************************* ! !! CSPDI sets the determinant and inverse of a complex symmetric packed matrix. ! ! ! Discussion: ! ! CSPDI uses the factors from cspfa. ! The matrix is stored in packed form. ! ! A division by zero will occur if the inverse is requested and CSPCO has ! set RCOND to 0.0 or CSPFA has set INFO nonzero. ! ! Parameters: ! ! on entry ! ! ap complex (n*(n+1)/2) ! the output from cspfa. ! ! Input, integer N, the order of the matrix. ! ! ipvt integer(n) ! the pivot vector from cspfa. ! ! work complex(n) ! work vector. contents ignored. ! ! job integer ! job has the decimal expansion ab where ! if b /= 0, the inverse is computed, ! if a /= 0, the determinant is computed, ! ! for example, job = 11 gives both. ! ! on return ! ! variables not requested by job are not used. ! ! 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. ! ! det complex(2) ! determinant of original matrix. ! determinant = det(1) * 10.0**det(2) ! with 1.0 <= abs(det(1)) < 10.0 ! or det(1) = 0.0. ! implicit none ! integer n ! complex ak complex akkp1 complex akp1 complex ap((n*(n+1))/2) real cabs1 complex cdotu complex d complex det(2) integer ij integer ik integer ikp1 integer iks 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 noinv complex t complex temp real, parameter :: ten = 10.0E+00 complex work(n) ! noinv = mod(job,10) == 0 nodet = mod(job,100)/10 == 0 if ( .not. nodet ) then det(1) = (1.0e0,0.0e0) det(2) = (0.0e0,0.0e0) t = (0.0e0,0.0e0) ik = 0 do k = 1, n kk = ik + k d = ap(kk) ! ! 2 by 2 block ! use det (d t) = (d/t * c - t) * t ! (t c) ! to avoid underflow/overflow troubles. ! take two passes through scaling. use t for flag. ! if (ipvt(k) <= 0) then if ( cabs1(t) == 0.0e0 ) then ikp1 = ik + k kkp1 = ikp1 + k t = ap(kkp1) d = (d/t)*ap(kkp1+1) - t else d = t t = (0.0e0,0.0e0) end if end if if ( .not. nodet ) then det(1) = d*det(1) if ( cabs1(det(1)) /= 0.0e0 ) then do while ( cabs1(det(1)) < 1.0e0 ) det(1) = cmplx(ten,0.0e0)*det(1) det(2) = det(2) - (1.0e0,0.0e0) end do do while (cabs1(det(1)) >= ten) det(1) = det(1)/cmplx(ten,0.0e0) det(2) = det(2) + (1.0e0,0.0e0) 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 if ( ipvt(k) >= 0) then ! ! 1 by 1 ! ap(kk) = (1.0e0,0.0e0)/ap(kk) if ( km1 >= 1 ) then call ccopy(km1,ap(ik+1),1,work,1) ij = 0 do j = 1, km1 jk = ik + j ap(jk) = cdotu(j,ap(ij+1),1,work,1) call caxpy(j-1,work(j),ap(ij+1),1,ap(ik+1),1) ij = ij + j end do ap(kk) = ap(kk) + cdotu(km1,work,1,ap(ik+1),1) end if kstep = 1 ! ! 2 by 2 ! else kkp1 = ikp1 + k t = ap(kkp1) ak = ap(kk)/t akp1 = ap(kkp1+1)/t akkp1 = ap(kkp1)/t d = t*(ak*akp1 - (1.0e0,0.0e0)) ap(kk) = akp1/d ap(kkp1+1) = ak/d ap(kkp1) = -akkp1/d if ( km1 >= 1) then call ccopy(km1,ap(ikp1+1),1,work,1) ij = 0 do j = 1, km1 jkp1 = ikp1 + j ap(jkp1) = cdotu(j,ap(ij+1),1,work,1) call caxpy(j-1,work(j),ap(ij+1),1,ap(ikp1+1),1) ij = ij + j end do ap(kkp1+1) = ap(kkp1+1) + cdotu(km1,work,1,ap(ikp1+1),1) ap(kkp1) = ap(kkp1) + cdotu(km1,ap(ik+1),1,ap(ikp1+1),1) call ccopy(km1,ap(ik+1),1,work,1) ij = 0 do j = 1, km1 jk = ik + j ap(jk) = cdotu(j,ap(ij+1),1,work,1) call caxpy(j-1,work(j),ap(ij+1),1,ap(ik+1),1) ij = ij + j end do ap(kk) = ap(kk) + cdotu(km1,work,1,ap(ik+1),1) 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),1,ap(ik+1),1) ksj = ik + ks do jb = ks, k j = k + ks - jb jk = ik + j temp = ap(jk) ap(jk) = 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) ik = ik + k + 1 k = k + kstep end do end if return end subroutine cspfa ( ap, n, ipvt, info ) ! !******************************************************************************* ! !! CSPFA factors a complex symmetric matrix stored in ! packed form by elimination with symmetric pivoting. ! ! to solve a*x = b , follow cspfa by cspsl. ! to compute inverse(a)*c , follow cspfa by cspsl. ! to compute determinant(a) , follow cspfa by cspdi. ! to compute inverse(a) , follow cspfa by cspdi. ! ! packed storage ! ! the following program segment will pack the upper ! triangle of a symmetric matrix. ! ! k = 0 ! do j = 1, n ! do i = 1, j ! k = k + 1 ! ap(k) = a(i,j) ! end do ! end do ! ! Parameters: ! ! on entry ! ! ap complex (n*(n+1)/2) ! the packed form of a symmetric matrix a . the ! columns of the upper triangle are stored sequentially ! in a one-dimensional array of length n*(n+1)/2 . ! see comments below for details. ! ! Input, integer N, the order of the matrix. ! ! on return ! ! ap 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*trans(u) ! where u is a product of permutation and unit ! upper triangular matrices , trans(u) is the ! transpose of u , and d is block diagonal ! with 1 by 1 and 2 by 2 blocks. ! ! ipvt integer(n) ! an integer vector of pivot indices. ! ! info integer ! = 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 cspsl or cspdi 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 imaxp1 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.0e0 + sqrt(17.0e0))/8.0e0 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.0e0) info = 1 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),1) imk = ik + imax colmax = cabs1(ap(imk)) if ( absakk >= alpha*colmax) then kstep = 1 swap = .false. ! ! Determine the largest off-diagonal element in row IMAX. ! else rowmax = 0.0e0 imaxp1 = imax + 1 im = imax*(imax - 1)/2 imj = im + 2*imax do j = imaxp1, k rowmax = max (rowmax,cabs1(ap(imj))) imj = imj + j end do if ( imax /= 1 ) then jmax = icamax(imax-1,ap(im+1),1) jmim = jmax + im rowmax = max (rowmax,cabs1(ap(jmim))) end if imim = imax + im if ( cabs1(ap(imim)) >= alpha*rowmax ) then kstep = 1 swap = .true. else if (absakk >= alpha*colmax*(colmax/rowmax)) 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.0e0) then ipvt(k) = k info = k ik = ik - (k - 1) if (kstep == 2) ik = ik - (k - 2) k = k - kstep cycle end if if ( kstep /= 2 ) then ! ! 1 x 1 pivot block. ! if ( swap ) then call cswap(imax,ap(im+1),1,ap(ik+1),1) imj = ik + imax do jj = imax, k j = k + imax - jj jk = ik + j t = ap(jk) ap(jk) = 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 = mulk call caxpy(j,t,ap(ik+1),1,ap(ij+1),1) ijj = ij + j ap(jk) = mulk ij = ij - (j - 1) end do ! ! Set the pivot array. ! ipvt(k) = k if (swap) ipvt(k) = imax ! ! 2 x 2 pivot block. ! else km1k = ik + k - 1 ikm1 = ik - (k - 1) if ( swap ) then call cswap(imax,ap(im+1),1,ap(ikm1+1),1) imj = ikm1 + imax do jj = imax, km1 j = km1 + imax - jj jkm1 = ikm1 + j t = ap(jkm1) ap(jkm1) = 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)/ap(km1k) denom = 1.0e0 - 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)/ap(km1k) mulk = (akm1*bk - bkm1)/denom mulkm1 = (ak*bkm1 - bk)/denom t = mulk call caxpy(j,t,ap(ik+1),1,ap(ij+1),1) t = mulkm1 call caxpy(j,t,ap(ikm1+1),1,ap(ij+1),1) ap(jk) = mulk ap(jkm1) = mulkm1 ijj = ij + j ij = ij - (j - 1) end do end if ! ! Set the pivot array. ! ipvt(k) = 1 - k if (swap) ipvt(k) = -imax ipvt(k-1) = ipvt(k) end if ik = ik - (k - 1) if (kstep == 2) ik = ik - (k - 2) k = k - kstep end do return end subroutine cspsl ( ap, n, ipvt, b ) ! !******************************************************************************* ! !! CSPSL solves the complex symmetric system ! a * x = b ! using the factors computed by cspfa. ! ! error condition ! ! a division by zero may occur if cspco has set rcond == 0.0 ! or cspfa has set info /= 0 . ! ! to compute inverse(a) * c where c is a matrix ! with p columns ! call cspfa(ap,n,ipvt,info) ! if (info /= 0) go to ... ! do j = 1, p ! call cspsl(ap,n,ipvt,c(1,j)) ! end do ! ! Parameters: ! ! on entry ! ! ap complex(n*(n+1)/2) ! the output from cspfa. ! ! Input, integer N, the order of the matrix. ! ! ipvt integer(n) ! the pivot vector from cspfa. ! ! 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 cdotu complex denom integer ik integer ikm1 integer ikp1 integer ipvt(n) integer k integer kk integer km1k integer km1km1 integer kp complex temp ! ! Loop backward applying the transformations and d inverse to b. ! k = n ik = (n*(n - 1))/2 do while ( k > 0 ) kk = ik + k if ( ipvt(k) >= 0 ) then ! ! 1 x 1 pivot block. ! if (k /= 1) then kp = ipvt(k) if ( kp /= k ) then call c_swap ( b(k), b(kp) ) end if call caxpy(k-1,b(k),ap(ik+1),1,b(1),1) end if ! ! Apply D inverse. ! b(k) = b(k)/ap(kk) k = k - 1 ik = ik - k ! ! 2 x 2 pivot block. ! else ikm1 = ik - (k - 1) if ( k /= 2 ) then kp = abs(ipvt(k)) if ( kp /= k - 1) then temp = b(k-1) b(k-1) = b(kp) b(kp) = temp end if call caxpy(k-2,b(k),ap(ik+1),1,b(1),1) call caxpy(k-2,b(k-1),ap(ikm1+1),1,b(1),1) end if ! ! Apply D inverse. ! km1k = ik + k - 1 kk = ik + k ak = ap(kk)/ap(km1k) km1km1 = ikm1 + k - 1 akm1 = ap(km1km1)/ap(km1k) bk = b(k)/ap(km1k) bkm1 = b(k-1)/ap(km1k) denom = ak*akm1 - 1.0e0 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 (ipvt(k) >= 0) then if ( k /= 1 ) then b(k) = b(k) + cdotu(k-1,ap(ik+1),1,b(1),1) kp = ipvt(k) if ( kp /= k ) then call c_swap ( b(k), b(kp) ) end if end if ik = ik + k k = k + 1 ! ! 2 x 2 pivot block. ! else if (k /= 1) then b(k) = b(k) + cdotu(k-1,ap(ik+1),1,b(1),1) ikp1 = ik + k b(k+1) = b(k+1) + cdotu(k-1,ap(ikp1+1),1,b(1),1) kp = abs(ipvt(k)) if (kp /= k) then temp = b(k) b(k) = b(kp) b(kp) = temp end if end if ik = ik + k + k + 1 k = k + 2 end if end do return end subroutine csvdc ( x, ldx, n, p, s, e, u, ldu, v, ldv, work, job, info ) ! !******************************************************************************* ! !! CSVDC is a subroutine to reduce a complex nxp matrix x by ! unitary transformations u and v to diagonal form. the ! diagonal elements s(i) are the singular values of x. the ! columns of u are the corresponding left singular vectors, ! and the columns of v the right singular vectors. ! ! Parameters: ! ! on entry ! ! x complex(ldx,p), where ldx>=n. ! x contains the matrix whose singular value ! decomposition is to be computed. x is ! destroyed by csvdc. ! ! Input, integer LDX, the leading dimension of X. ! ! Input, integer N, the number of rows of the matrix. ! ! p integer. ! p is the number of columns of the matrix x. ! ! Input, integer LDU, the leading dimension of U. ! ! Input, integer LDV, the leading dimension of V. ! ! work complex(n). ! work is a scratch array. ! ! job integer. ! job controls the computation of the singular ! vectors. it has the decimal expansion ab ! with the following meaning ! ! a==0 do not compute the left singular ! vectors. ! a==1 return the n left singular vectors ! in u. ! a>=2 returns the first min(n,p) ! left singular vectors in u. ! b==0 do not compute the right singular ! vectors. ! b==1 return the right singular vectors ! in v. ! ! on return ! ! s complex(mm), where mm=min(n+1,p). ! the first min(n,p) entries of s contain the ! singular values of x arranged in descending ! order of magnitude. ! ! e complex(p). ! e ordinarily contains zeros. however see the ! discussion of info for exceptions. ! ! u complex(ldu,k), where ldu>=n. if joba==1 then ! k==n, if joba>=2 then ! k==min(n,p). ! u contains the matrix of left singular vectors. ! u is not referenced if joba==0. if n<=p ! or if joba>2, then u may be identified with x ! in the subroutine call. ! ! v complex(ldv,p), where ldv>=p. ! v contains the matrix of right singular vectors. ! v is not referenced if jobb==0. if p<=n, ! then v may be identified whth x in the ! subroutine call. ! ! info integer. ! the singular values (and their corresponding ! singular vectors) s(info+1),s(info+2),...,s(m) ! are correct (here m=min(n,p)). thus if ! info==0, all the singular values and their ! vectors are correct. in any event, the matrix ! b = hermitian(u)*x*v is the bidiagonal matrix ! with the elements of s on its diagonal and the ! elements of e on its super-diagonal (hermitian(u) ! is the conjugate-transpose of u). thus the ! singular values of x and b are the same. ! implicit none ! integer ldu integer ldv integer ldx integer n integer p ! real b real c real cabs1 complex cdotc real cs complex csign complex e(p) real el real emm1 real f real g integer i integer info integer iter integer j integer job integer jobu integer k integer kase integer kk integer l integer ll integer lls integer lm1 integer lp1 integer ls integer lu integer m integer, parameter :: maxit = 30 integer mm integer mm1 integer mp1 integer nct integer nctp1 integer ncu integer nrt integer nrtp1 complex r complex s(*) real scnrm2 real scale real shift real sl real sm real smm1 real sn complex t real t1 real test complex u(ldu,*) complex v(ldv,p) logical wantu logical wantv complex work(n) complex x(ldx,p) real ztest ! ! Determine what is to be computed. ! wantu = .false. wantv = .false. jobu = mod(job,100)/10 ncu = n if (jobu > 1) ncu = min(n,p) if (jobu /= 0) wantu = .true. if (mod(job,10) /= 0) wantv = .true. ! ! Reduce X to bidiagonal form, storing the diagonal elements ! in S and the super-diagonal elements in E. ! info = 0 nct = min(n-1,p) nrt = max(0,min(p-2,n)) lu = max(nct,nrt) do l = 1, lu lp1 = l + 1 ! ! Compute the transformation for the L-th column and ! place the L-th diagonal in S(L). ! if (l <= nct) then s(l) = cmplx(scnrm2(n-l+1,x(l,l),1),0.0e0) if (cabs1(s(l)) /= 0.0e0) then if (cabs1(x(l,l)) /= 0.0e0) s(l) = csign(s(l),x(l,l)) call cscal(n-l+1,1.0e0/s(l),x(l,l),1) x(l,l) = (1.0e0,0.0e0) + x(l,l) end if s(l) = -s(l) end if do j = lp1, p if (l <= nct) then if (cabs1(s(l)) /= 0.0e0) then t = -cdotc(n-l+1,x(l,l),1,x(l,j),1)/x(l,l) call caxpy(n-l+1,t,x(l,l),1,x(l,j),1) end if end if ! ! Place the L-th row of X into E for the ! subsequent calculation of the row transformation. ! e(j) = conjg(x(l,j)) end do ! ! Place the transformation in U for subsequent back multiplication. ! if ( wantu .and. l <= nct ) then u(l:n,l) = x(l:n,l) end if if ( l <= nrt ) then ! ! Compute the L-th row transformation and place the ! L-th super-diagonal in E(L). ! e(l) = cmplx(scnrm2(p-l,e(lp1),1),0.0e0) if ( cabs1(e(l)) /= 0.0e0 ) then if (cabs1(e(lp1)) /= 0.0e0) e(l) = csign(e(l),e(lp1)) call cscal(p-l,1.0e0/e(l),e(lp1),1) e(lp1) = (1.0e0,0.0e0) + e(lp1) end if e(l) = -conjg(e(l)) ! ! Apply the transformation. ! if (lp1 <= n .and. cabs1(e(l)) /= 0.0e0) then work(lp1:n) = (0.0e0,0.0e0) do j = lp1, p call caxpy(n-l,e(j),x(lp1,j),1,work(lp1),1) end do do j = lp1, p call caxpy(n-l,conjg(-e(j)/e(lp1)),work(lp1),1, & x(lp1,j),1) end do end if ! ! Place the transformation in V for subsequent back multiplication. ! if ( wantv ) then do i = lp1, p v(i,l) = e(i) end do end if end if end do ! ! Set up the final bidiagonal matrix of order M. ! m = min(p,n+1) nctp1 = nct + 1 nrtp1 = nrt + 1 if (nct < p) s(nctp1) = x(nctp1,nctp1) if (n < m) s(m) = (0.0e0,0.0e0) if (nrtp1 < m) e(nrtp1) = x(nrtp1,m) e(m) = (0.0e0,0.0e0) ! ! If required, generate U. ! if ( wantu ) then do j = nctp1, ncu u(1:n,j) = (0.0e0,0.0e0) u(j,j) = (1.0e0,0.0e0) end do do ll = 1, nct l = nct - ll + 1 if (cabs1(s(l)) /= 0.0e0) then lp1 = l + 1 do j = l+1, ncu t = -cdotc(n-l+1,u(l,l),1,u(l,j),1)/u(l,l) call caxpy(n-l+1,t,u(l,l),1,u(l,j),1) end do call cscal(n-l+1,(-1.0e0,0.0e0),u(l,l),1) u(l,l) = (1.0e0,0.0e0) + u(l,l) u(1:l-1,l) = (0.0e0,0.0e0) else u(1:n,l) = (0.0e0,0.0e0) u(l,l) = (1.0e0,0.0e0) end if end do end if ! ! If it is required, generate V. ! if ( wantv ) then do ll = 1, p l = p - ll + 1 lp1 = l + 1 if (l <= nrt) then if ( cabs1(e(l)) /= 0.0e0) then do j = lp1, p t = -cdotc(p-l,v(lp1,l),1,v(lp1,j),1)/v(lp1,l) call caxpy(p-l,t,v(lp1,l),1,v(lp1,j),1) end do end if end if v(1:p,l) = (0.0e0,0.0e0) v(l,l) = (1.0e0,0.0e0) end do end if ! ! Transform S and E so that they are real. ! do i = 1, m if ( cabs1(s(i)) /= 0.0e0 ) then t = cmplx(cabs(s(i)),0.0e0) r = s(i)/t s(i) = t if (i < m) e(i) = e(i)/r if (wantu) call cscal(n,r,u(1,i),1) end if if ( i == m ) then exit end if if ( cabs1(e(i)) /= 0.0e0 ) then t = cmplx(cabs(e(i)),0.0e0) r = t/e(i) e(i) = t s(i+1) = s(i+1)*r if (wantv) call cscal(p,r,v(1,i+1),1) end if end do ! ! Main iteration loop for the singular values. ! mm = m iter = 0 do ! ! Quit if all the singular values have been found. ! if ( m == 0 ) then exit end if ! ! If too many iterations have been performed, set flag and return. ! if ( iter >= maxit ) then info = m exit end if ! ! This section of the program inspects for negligible elements in S and E. ! ! On completion the variables kase and l are set as follows. ! ! kase = 1 if s(m) and e(l-1) are negligible and l= real(s(l+1)) ) then exit end if t = s(l) s(l) = s(l+1) s(l+1) = t if ( wantv .and. l < p ) then call cswap(p,v(1,l),1,v(1,l+1),1) end if if ( wantu .and. l < n ) then call cswap(n,u(1,l),1,u(1,l+1),1) end if l = l + 1 end do iter = 0 m = m - 1 end if end do return end subroutine ctrco ( t, ldt, n, rcond, z, job ) ! !******************************************************************************* ! !! CTRCO estimates the condition of a complex triangular matrix. ! ! ! Parameters: ! ! on entry ! ! t complex(ldt,n) ! t contains the triangular matrix. the zero ! elements of the matrix are not referenced, and ! the corresponding elements of the array can be ! used to store other information. ! ! Input, integer LDT, the leading dimension of T. ! ! Input, integer N, the order of the matrix. ! ! Input, integer JOB, indicates if matrix is upper or lower triangular. ! 0, lower triangular. ! nonzero, upper triangular. ! ! on return ! ! rcond real ! an estimate of the reciprocal condition of t . ! for the system t*x = b , relative perturbations ! in t 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 t may be singular to working ! precision. in particular, rcond is zero if ! exact singularity is detected or the estimate ! underflows. ! ! z complex(n) ! a work vector whose contents are usually unimportant. ! if t 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 ldt integer n ! real cabs1 complex csign1 complex ek integer i1 integer j integer j1 integer j2 integer job integer k integer kk integer l logical lower real rcond real s real scasum real sm complex t(ldt,n) real tnorm complex w complex wk complex wkm real ynorm complex z(n) ! lower = job == 0 ! ! Compute 1-norm of T ! tnorm = 0.0e0 do j = 1, n l = j if (lower) l = n + 1 - j i1 = 1 if (lower) i1 = j tnorm = max (tnorm,scasum(l,t(i1,j),1)) end do ! ! RCOND = 1/(norm(t)*(estimate of norm(inverse(t)))). ! ! Estimate = norm(z)/norm(y) where t*z = y and hermitian(t)*y = E. ! ! hermitian(t) is the conjugate transpose of T. ! ! The components of E are chosen to cause maximum local ! growth in the elements of Y. ! ! The vectors are frequently rescaled to avoid overflow. ! ! Solve hermitian(t)*y = E. ! ek = (1.0e0,0.0e0) z(1:n) = (0.0e0,0.0e0) do kk = 1, n k = kk if (lower) k = n + 1 - kk if (cabs1(z(k)) /= 0.0e0) ek = csign1(ek,-z(k)) if ( cabs1(ek-z(k)) > cabs1(t(k,k)) ) then s = cabs1(t(k,k))/cabs1(ek-z(k)) call csscal(n,s,z,1) ek = cmplx(s,0.0e0)*ek end if wk = ek - z(k) wkm = -ek - z(k) s = cabs1(wk) sm = cabs1(wkm) if ( cabs1(t(k,k)) /= 0.0e0 ) then wk = wk/conjg(t(k,k)) wkm = wkm/conjg(t(k,k)) else wk = (1.0e0,0.0e0) wkm = (1.0e0,0.0e0) end if if ( kk /= n ) then j1 = k + 1 if (lower) j1 = 1 j2 = n if (lower) j2 = k - 1 do j = j1, j2 sm = sm + cabs1(z(j)+wkm*conjg(t(k,j))) z(j) = z(j) + wk*conjg(t(k,j)) s = s + cabs1(z(j)) end do if ( s < sm ) then w = wkm - wk wk = wkm do j = j1, j2 z(j) = z(j) + w*conjg(t(k,j)) end do end if end if z(k) = wk end do s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = 1.0e0 ! ! Solve t*z = Y. ! do kk = 1, n k = n + 1 - kk if (lower) k = kk if ( cabs1(z(k)) > cabs1(t(k,k)) ) then s = cabs1(t(k,k))/cabs1(z(k)) call csscal(n,s,z,1) ynorm = s*ynorm end if if (cabs1(t(k,k)) /= 0.0e0) z(k) = z(k)/t(k,k) if (cabs1(t(k,k)) == 0.0e0) z(k) = (1.0e0,0.0e0) i1 = 1 if (lower) i1 = k + 1 if ( kk < n ) then w = -z(k) call caxpy(n-kk,w,t(i1,k),1,z(i1),1) end if end do ! ! Make ZNORM = 1. ! s = 1.0e0/scasum(n,z,1) call csscal(n,s,z,1) ynorm = s*ynorm if ( tnorm /= 0.0e0 ) then rcond = ynorm / tnorm else rcond = 0.0e0 end if return end subroutine ctrdi ( t, ldt, n, det, job, info ) ! !******************************************************************************* ! !! CTRDI computes the determinant and inverse of a complex triangular matrix. ! ! ! Parameters: ! ! on entry ! ! t complex(ldt,n) ! t contains the triangular matrix. the zero ! elements of the matrix are not referenced, and ! the corresponding elements of the array can be ! used to store other information. ! ! Input, integer LDT, the leading dimension of T. ! ! Input, integer N, the order of the matrix. ! ! job integer ! = 010 no det, inverse of lower triangular. ! = 011 no det, inverse of upper triangular. ! = 100 det, no inverse. ! = 110 det, inverse of lower triangular. ! = 111 det, inverse of upper triangular. ! ! on return ! ! t inverse of original matrix if requested. ! otherwise unchanged. ! ! det complex(2) ! determinant of 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 . ! ! info integer ! info contains zero if the system is nonsingular ! and the inverse is requested. ! otherwise info contains the index of ! a zero diagonal element of t. ! ! implicit none ! integer ldt integer n ! real cabs1 complex det(2) integer i integer info integer j integer job integer k integer kb integer km1 integer kp1 complex t(ldt,n) complex temp real, parameter :: ten = 10.0E+00 ! if ( job/100 /= 0 ) then det(1) = (1.0e0,0.0e0) det(2) = (0.0e0,0.0e0) do i = 1, n det(1) = t(i,i)*det(1) if ( cabs1(det(1)) == 0.0e0 ) then exit end if do while ( cabs1(det(1)) < 1.0e0 ) det(1) = cmplx(ten,0.0e0)*det(1) det(2) = det(2) - (1.0e0,0.0e0) end do do while ( cabs1(det(1)) >= ten ) det(1) = det(1)/cmplx(ten,0.0e0) det(2) = det(2) + (1.0e0,0.0e0) end do end do end if ! ! Compute inverse of upper triangular matrix. ! if ( mod ( job/10, 10 ) /= 0 ) then if (mod(job,10) /= 0) then info = 0 do k = 1, n if ( cabs1(t(k,k)) == 0.0e0 ) then info = k exit end if t(k,k) = (1.0e0,0.0e0)/t(k,k) temp = -t(k,k) call cscal(k-1,temp,t(1,k),1) do j = k+1, n temp = t(k,j) t(k,j) = (0.0e0,0.0e0) call caxpy(k,temp,t(1,k),1,t(1,j),1) end do end do ! ! Compute inverse of lower triangular matrix. ! else info = 0 do kb = 1, n k = n + 1 - kb if ( cabs1(t(k,k)) == 0.0e0 ) then info = k exit end if t(k,k) = (1.0e0,0.0e0)/t(k,k) if (k /= n) then temp = -t(k,k) call cscal(n-k,temp,t(k+1,k),1) end if do j = 1, k-1 temp = t(k,j) t(k,j) = (0.0e0,0.0e0) call caxpy(n-k+1,temp,t(k,k),1,t(k,j),1) end do end do end if end if return end subroutine ctrsl ( t, ldt, n, b, job, info ) ! !******************************************************************************* ! !! CTRSL solves triangular systems ! ! of the form ! ! t * x = b ! or ! hermitian(t) * x = b ! ! where t is a triangular matrix of order n. here hermitian(t) ! denotes the conjugate transpose of the matrix t. ! ! Parameters: ! ! on entry ! ! t complex(ldt,n) ! t contains the matrix of the system. the zero ! elements of the matrix are not referenced, and ! the corresponding elements of the array can be ! used to store other information. ! ! Input, integer LDT, the leading dimension of T. ! ! Input, integer N, the order of the matrix. ! ! Input/output, complex B(N). On input, the right hand side. ! On output, the solution. ! ! job integer ! job specifies what kind of system is to be solved. ! if job is ! ! 00 solve t*x=b, t lower triangular, ! 01 solve t*x=b, t upper triangular, ! 10 solve hermitian(t)*x=b, t lower triangular, ! 11 solve hermitian(t)*x=b, t upper triangular. ! ! on return ! ! info integer ! info contains zero if the system is nonsingular. ! otherwise info contains the index of ! the first zero diagonal element of t. ! implicit none ! integer ldt integer n ! complex b(n) real cabs1 integer case complex cdotc integer info integer j integer jj integer job complex t(ldt,n) complex temp ! ! Check for zero diagonal elements. ! do info = 1, n if ( cabs1 ( t(info,info) ) == 0.0e0 ) then return end if end do info = 0 ! ! Determine the task and go to it. ! case = 1 if (mod(job,10) /= 0) case = 2 if (mod(job,100)/10 /= 0) case = case + 2 ! ! Solve t*x=b for T lower triangular. ! if ( case == 1 ) then b(1) = b(1)/t(1,1) do j = 2, n temp = -b(j-1) call caxpy(n-j+1,temp,t(j,j-1),1,b(j),1) b(j) = b(j)/t(j,j) end do ! ! Solve t*x=b for T upper triangular. ! else if ( case == 2 ) then b(n) = b(n)/t(n,n) do jj = 2, n j = n - jj + 1 temp = -b(j+1) call caxpy(j,temp,t(1,j+1),1,b(1),1) b(j) = b(j)/t(j,j) end do ! ! Solve hermitian(t)*x=b for T lower triangular. ! else if ( case == 3 ) then b(n) = b(n)/conjg(t(n,n)) do jj = 2, n j = n - jj + 1 b(j) = b(j) - cdotc(jj-1,t(j+1,j),1,b(j+1),1) b(j) = b(j)/conjg(t(j,j)) end do ! ! Solve hermitian(t)*x=b for T upper triangular. ! else if ( case == 4 ) then b(1) = b(1)/conjg(t(1,1)) do j = 2, n b(j) = b(j) - cdotc(j-1,t(1,j),1,b(1),1) b(j) = b(j)/conjg(t(j,j)) end do end if return end subroutine timestamp ( ) ! !******************************************************************************* ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! ! Example: ! ! May 31 2001 9:45:54.872 AM ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none ! character ( len = 8 ) ampm integer d character ( len = 8 ) date integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s character ( len = 10 ) time integer values(8) integer y character ( len = 5 ) zone ! call date_and_time ( date, time, zone, values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end