subroutine schdc(a,lda,p,work,jpvt,job,info) c*********************************************************************72 c c schdc computes the cholesky decomposition of a positive definite c matrix. a pivoting option allows the user to estimate the c condition of a positive definite matrix or determine the rank c of a positive semidefinite matrix. c c on entry c c a real(lda,p). c a contains the matrix whose decomposition is to c be computed. onlt the upper half of a need be stored. c the lower part of the array a is not referenced. c c lda integer. c lda is the leading dimension of the array a. c c p integer. c p is the order of the matrix. c c work real. c work is a work array. c c jpvt integer(p). c jpvt contains integers that control the selection c of the pivot elements, if pivoting has been requested. c each diagonal element a(k,k) c is placed in one of three classes according to the c value of jpvt(k). c c if jpvt(k) .gt. 0, then x(k) is an initial c element. c c if jpvt(k) .eq. 0, then x(k) is a free element. c c if jpvt(k) .lt. 0, then x(k) is a final element. c c before the decomposition is computed, initial elements c are moved by symmetric row and column interchanges to c the beginning of the array a and final c elements to the end. both initial and final elements c are frozen in place during the computation and only c free elements are moved. at the k-th stage of the c reduction, if a(k,k) is occupied by a free element c it is interchanged with the largest free element c a(l,l) with l .ge. k. jpvt is not referenced if c job .eq. 0. c c job integer. c job is an integer that initiates column pivoting. c if job .eq. 0, no pivoting is done. c if job .ne. 0, pivoting is done. c c on return c c a a contains in its upper half the cholesky factor c of the matrix a as it has been permuted by pivoting. c c jpvt jpvt(j) contains the index of the diagonal element c of a that was moved into the j-th position, c provided pivoting was requested. c c info contains the index of the last positive diagonal c element of the cholesky factor. c c for positive definite matrices info = p is the normal return. c for pivoting with positive semidefinite matrices info will c in general be less than p. however, info may be greater than c the rank of a, since rounding error can cause an otherwise zero c element to be positive. indefinite systems will always cause c info to be less than p. c c linpack. this version dated 03/19/79 . c j.j. dongarra and g.w. stewart, argonne national laboratory and c university of maryland. c c c blas saxpy,sswap c fortran sqrt c integer lda,p,jpvt(*),job,info real a(lda,1),work(*) integer pu,pl,plp1,j,jp,jt,k,kb,km1,kp1,l,maxl real temp real maxdia logical swapk,negk c pl = 1 pu = 0 info = p if (job .eq. 0) go to 160 c c pivoting has been requested. rearrange the c the elements according to jpvt. c do 70 k = 1, p swapk = jpvt(k) .gt. 0 negk = jpvt(k) .lt. 0 jpvt(k) = k if (negk) jpvt(k) = -jpvt(k) if (.not.swapk) go to 60 if (k .eq. pl) go to 50 call sswap(pl-1,a(1,k),1,a(1,pl),1) temp = a(k,k) a(k,k) = a(pl,pl) a(pl,pl) = temp plp1 = pl + 1 if (p .lt. plp1) go to 40 do 30 j = plp1, p if (j .ge. k) go to 10 temp = a(pl,j) a(pl,j) = a(j,k) a(j,k) = temp go to 20 10 continue if (j .eq. k) go to 20 temp = a(k,j) a(k,j) = a(pl,j) a(pl,j) = temp 20 continue 30 continue 40 continue jpvt(k) = jpvt(pl) jpvt(pl) = k 50 continue pl = pl + 1 60 continue 70 continue pu = p if (p .lt. pl) go to 150 do 140 kb = pl, p k = p - kb + pl if (jpvt(k) .ge. 0) go to 130 jpvt(k) = -jpvt(k) if (pu .eq. k) go to 120 call sswap(k-1,a(1,k),1,a(1,pu),1) temp = a(k,k) a(k,k) = a(pu,pu) a(pu,pu) = temp kp1 = k + 1 if (p .lt. kp1) go to 110 do 100 j = kp1, p if (j .ge. pu) go to 80 temp = a(k,j) a(k,j) = a(j,pu) a(j,pu) = temp go to 90 80 continue if (j .eq. pu) go to 90 temp = a(k,j) a(k,j) = a(pu,j) a(pu,j) = temp 90 continue 100 continue 110 continue jt = jpvt(k) jpvt(k) = jpvt(pu) jpvt(pu) = jt 120 continue pu = pu - 1 130 continue 140 continue 150 continue 160 continue do 270 k = 1, p c c reduction loop. c maxdia = a(k,k) kp1 = k + 1 maxl = k c c determine the pivot element. c if (k .lt. pl .or. k .ge. pu) go to 190 do 180 l = kp1, pu if (a(l,l) .le. maxdia) go to 170 maxdia = a(l,l) maxl = l 170 continue 180 continue 190 continue c c quit if the pivot element is not positive. c if (maxdia .gt. 0.0e0) go to 200 info = k - 1 c ......exit go to 280 200 continue if (k .eq. maxl) go to 210 c c start the pivoting and update jpvt. c km1 = k - 1 call sswap(km1,a(1,k),1,a(1,maxl),1) a(maxl,maxl) = a(k,k) a(k,k) = maxdia jp = jpvt(maxl) jpvt(maxl) = jpvt(k) jpvt(k) = jp 210 continue c c reduction step. pivoting is contained across the rows. c work(k) = sqrt(a(k,k)) a(k,k) = work(k) if (p .lt. kp1) go to 260 do 250 j = kp1, p if (k .eq. maxl) go to 240 if (j .ge. maxl) go to 220 temp = a(k,j) a(k,j) = a(j,maxl) a(j,maxl) = temp go to 230 220 continue if (j .eq. maxl) go to 230 temp = a(k,j) a(k,j) = a(maxl,j) a(maxl,j) = temp 230 continue 240 continue a(k,j) = a(k,j)/work(k) work(j) = a(k,j) temp = -a(k,j) call saxpy(j-k,temp,work(kp1),1,a(kp1,j),1) 250 continue 260 continue 270 continue 280 continue return end subroutine schdd(r,ldr,p,x,z,ldz,nz,y,rho,c,s,info) c*********************************************************************72 integer ldr,p,ldz,nz,info real r(ldr,1),x(1),z(ldz,1),y(1),s(1) real rho(1),c(1) c c schdd downdates an augmented cholesky decomposition or the c triangular factor of an augmented qr decomposition. c specifically, given an upper triangular matrix r of order p, a c row vector x, a column vector z, and a scalar y, schdd c determineds a orthogonal matrix u and a scalar zeta such that c c (r z ) (rr zz) c u * ( ) = ( ) , c (0 zeta) ( x y) c c where rr is upper triangular. if r and z have been obtained c from the factorization of a least squares problem, then c rr and zz are the factors corresponding to the problem c with the observation (x,y) removed. in this case, if rho c is the norm of the residual vector, then the norm of c the residual vector of the downdated problem is c sqrt(rho**2 - zeta**2). schdd will simultaneously downdate c several triplets (z,y,rho) along with r. c for a less terse description of what schdd does and how c it may be applied, see the linpack guide. c c the matrix u is determined as the product u(1)*...*u(p) c where u(i) is a rotation in the (p+1,i)-plane of the c form c c ( c(i) -s(i) ) c ( ) . c ( s(i) c(i) ) c c the rotations are chosen so that c(i) is real. c c the user is warned that a given downdating problem may c be impossible to accomplish or may produce c inaccurate results. for example, this can happen c if x is near a vector whose removal will reduce the c rank of r. beware. c c on entry c c r real(ldr,p), where ldr .ge. p. c r contains the upper triangular matrix c that is to be downdated. the part of r c below the diagonal is not referenced. c c ldr integer. c ldr is the leading dimension fo the array r. c c p integer. c p is the order of the matrix r. c c x real(p). c x contains the row vector that is to c be removed from r. x is not altered by schdd. c c z real(ldz,nz), where ldz .ge. p. c z is an array of nz p-vectors which c are to be downdated along with r. c c ldz integer. c ldz is the leading dimension of the array z. c c nz integer. c nz is the number of vectors to be downdated c nz may be zero, in which case z, y, and rho c are not referenced. c c y real(nz). c y contains the scalars for the downdating c of the vectors z. y is not altered by schdd. c c rho real(nz). c rho contains the norms of the residual c vectors that are to be downdated. c c on return c c r c z contain the downdated quantities. c rho c c c real(p). c c contains the cosines of the transforming c rotations. c c s real(p). c s contains the sines of the transforming c rotations. c c info integer. c info is set as follows. c c info = 0 if the entire downdating c was successful. c c info =-1 if r could not be downdated. c in this case, all quantities c are left unaltered. c c info = 1 if some rho could not be c downdated. the offending rhos are c set to -1. c c linpack. this version dated 08/14/78 . c g.w. stewart, university of maryland, argonne national lab. c c schdd uses the following functions and subprograms. c c fortran abs c blas sdot, snrm2 c integer i,ii,j real a,alpha,azeta,norm,snrm2 real sdot,t,zeta,b,xx c c solve the system trans(r)*a = x, placing the result c in the array s. c info = 0 s(1) = x(1)/r(1,1) if (p .lt. 2) go to 20 do 10 j = 2, p s(j) = x(j) - sdot(j-1,r(1,j),1,s,1) s(j) = s(j)/r(j,j) 10 continue 20 continue norm = snrm2(p,s,1) if (norm .lt. 1.0e0) go to 30 info = -1 go to 120 30 continue alpha = sqrt(1.0e0-norm**2) c c determine the transformations. c do 40 ii = 1, p i = p - ii + 1 scale = alpha + abs(s(i)) a = alpha/scale b = s(i)/scale norm = sqrt(a**2+b**2+0.0e0**2) c(i) = a/norm s(i) = b/norm alpha = scale*norm 40 continue c c apply the transformations to r. c do 60 j = 1, p xx = 0.0e0 do 50 ii = 1, j i = j - ii + 1 t = c(i)*xx + s(i)*r(i,j) r(i,j) = c(i)*r(i,j) - s(i)*xx xx = t 50 continue 60 continue c c if required, downdate z and rho. c if (nz .lt. 1) go to 110 do 100 j = 1, nz zeta = y(j) do 70 i = 1, p z(i,j) = (z(i,j) - s(i)*zeta)/c(i) zeta = c(i)*zeta - s(i)*z(i,j) 70 continue azeta = abs(zeta) if (azeta .le. rho(j)) go to 80 info = 1 rho(j) = -1.0e0 go to 90 80 continue rho(j) = rho(j)*sqrt(1.0e0-(azeta/rho(j))**2) 90 continue 100 continue 110 continue 120 continue return end subroutine schex(r,ldr,p,k,l,z,ldz,nz,c,s,job) c*********************************************************************72 integer ldr,p,k,l,ldz,nz,job real r(ldr,1),z(ldz,1),s(1) real c(1) c c schex updates the cholesky factorization c c a = trans(r)*r c c of a positive definite matrix a of order p under diagonal c permutations of the form c c trans(e)*a*e c c where e is a permutation matrix. specifically, given c an upper triangular matrix r and a permutation matrix c e (which is specified by k, l, and job), schex determines c a orthogonal matrix u such that c c u*r*e = rr, c c where rr is upper triangular. at the users option, the c transformation u will be multiplied into the array z. c if a = trans(x)*x, so that r is the triangular part of the c qr factorization of x, then rr is the triangular part of the c qr factorization of x*e, i.e. x with its columns permuted. c for a less terse description of what schex does and how c it may be applied, see the linpack guide. c c the matrix q is determined as the product u(l-k)*...*u(1) c of plane rotations of the form c c ( c(i) s(i) ) c ( ) , c ( -s(i) c(i) ) c c where c(i) is real, the rows these rotations operate on c are described below. c c there are two types of permutations, which are determined c by the value of job. c c 1. right circular shift (job = 1). c c the columns are rearranged in the following order. c c 1,...,k-1,l,k,k+1,...,l-1,l+1,...,p. c c u is the product of l-k rotations u(i), where u(i) c acts in the (l-i,l-i+1)-plane. c c 2. left circular shift (job = 2). c the columns are rearranged in the following order c c 1,...,k-1,k+1,k+2,...,l,k,l+1,...,p. c c u is the product of l-k rotations u(i), where u(i) c acts in the (k+i-1,k+i)-plane. c c on entry c c r real(ldr,p), where ldr.ge.p. c r contains the upper triangular factor c that is to be updated. elements of r c below the diagonal are not referenced. c c ldr integer. c ldr is the leading dimension of the array r. c c p integer. c p is the order of the matrix r. c c k integer. c k is the first column to be permuted. c c l integer. c l is the last column to be permuted. c l must be strictly greater than k. c c z real(ldz,nz), where ldz.ge.p. c z is an array of nz p-vectors into which the c transformation u is multiplied. z is c not referenced if nz = 0. c c ldz integer. c ldz is the leading dimension of the array z. c c nz integer. c nz is the number of columns of the matrix z. c c job integer. c job determines the type of permutation. c job = 1 right circular shift. c job = 2 left circular shift. c c on return c c r contains the updated factor. c c z contains the updated matrix z. c c c real(p). c c contains the cosines of the transforming rotations. c c s real(p). c s contains the sines of the transforming rotations. c c linpack. this version dated 08/14/78 . c g.w. stewart, university of maryland, argonne national lab. c c schex uses the following functions and subroutines. c c blas srotg c fortran min0 c integer i,ii,il,iu,j,jj,km1,kp1,lmk,lm1 real t c c initialize c km1 = k - 1 kp1 = k + 1 lmk = l - k lm1 = l - 1 c c perform the appropriate task. c go to (10,130), job c c right circular shift. c 10 continue c c reorder the columns. c do 20 i = 1, l ii = l - i + 1 s(i) = r(ii,l) 20 continue do 40 jj = k, lm1 j = lm1 - jj + k do 30 i = 1, j r(i,j+1) = r(i,j) 30 continue r(j+1,j+1) = 0.0e0 40 continue if (k .eq. 1) go to 60 do 50 i = 1, km1 ii = l - i + 1 r(i,k) = s(ii) 50 continue 60 continue c c calculate the rotations. c t = s(1) do 70 i = 1, lmk call srotg(s(i+1),t,c(i),s(i)) t = s(i+1) 70 continue r(k,k) = t do 90 j = kp1, p il = max0(1,l-j+1) do 80 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) - s(ii)*r(i,j) r(i,j) = t 80 continue 90 continue c c if required, apply the transformations to z. c if (nz .lt. 1) go to 120 do 110 j = 1, nz do 100 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) - s(ii)*z(i,j) z(i,j) = t 100 continue 110 continue 120 continue go to 260 c c left circular shift c 130 continue c c reorder the columns c do 140 i = 1, k ii = lmk + i s(ii) = r(i,k) 140 continue do 160 j = k, lm1 do 150 i = 1, j r(i,j) = r(i,j+1) 150 continue jj = j - km1 s(jj) = r(j+1,j+1) 160 continue do 170 i = 1, k ii = lmk + i r(i,l) = s(ii) 170 continue do 180 i = kp1, l r(i,l) = 0.0e0 180 continue c c reduction loop. c do 220 j = k, p if (j .eq. k) go to 200 c c apply the rotations. c iu = min0(j-1,l-1) do 190 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) - s(ii)*r(i,j) r(i,j) = t 190 continue 200 continue if (j .ge. l) go to 210 jj = j - k + 1 t = s(jj) call srotg(r(j,j),t,c(jj),s(jj)) 210 continue 220 continue c c apply the rotations to z. c if (nz .lt. 1) go to 250 do 240 j = 1, nz do 230 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) - s(ii)*z(i,j) z(i,j) = t 230 continue 240 continue 250 continue 260 continue return end subroutine schud(r,ldr,p,x,z,ldz,nz,y,rho,c,s) c*********************************************************************72 integer ldr,p,ldz,nz real rho(1),c(1) real r(ldr,1),x(1),z(ldz,1),y(1),s(1) c c schud updates an augmented cholesky decomposition of the c triangular part of an augmented qr decomposition. specifically, c given an upper triangular matrix r of order p, a row vector c x, a column vector z, and a scalar y, schud determines a c untiary matrix u and a scalar zeta such that c c c (r z) (rr zz ) c u * ( ) = ( ) , c (x y) ( 0 zeta) c c where rr is upper triangular. if r and z have been c obtained from the factorization of a least squares c problem, then rr and zz are the factors corresponding to c the problem with the observation (x,y) appended. in this c case, if rho is the norm of the residual vector, then the c norm of the residual vector of the updated problem is c sqrt(rho**2 + zeta**2). schud will simultaneously update c several triplets (z,y,rho). c for a less terse description of what schud does and how c it may be applied, see the linpack guide. c c the matrix u is determined as the product u(p)*...*u(1), c where u(i) is a rotation in the (i,p+1) plane of the c form c c ( c(i) s(i) ) c ( ) . c ( -s(i) c(i) ) c c the rotations are chosen so that c(i) is real. c c on entry c c r real(ldr,p), where ldr .ge. p. c r contains the upper triangular matrix c that is to be updated. the part of r c below the diagonal is not referenced. c c ldr integer. c ldr is the leading dimension of the array r. c c p integer. c p is the order of the matrix r. c c x real(p). c x contains the row to be added to r. x is c not altered by schud. c c z real(ldz,nz), where ldz .ge. p. c z is an array containing nz p-vectors to c be updated with r. c c ldz integer. c ldz is the leading dimension of the array z. c c nz integer. c nz is the number of vectors to be updated c nz may be zero, in which case z, y, and rho c are not referenced. c c y real(nz). c y contains the scalars for updating the vectors c z. y is not altered by schud. c c rho real(nz). c rho contains the norms of the residual c vectors that are to be updated. if rho(j) c is negative, it is left unaltered. c c on return c c rc c rho contain the updated quantities. c z c c c real(p). c c contains the cosines of the transforming c rotations. c c s real(p). c s contains the sines of the transforming c rotations. c c linpack. this version dated 08/14/78 . c g.w. stewart, university of maryland, argonne national lab. c c schud uses the following functions and subroutines. c c extended blas srotg c fortran sqrt c integer i,j,jm1 real azeta,scale real t,xj,zeta c c update r. c do 30 j = 1, p xj = x(j) c c apply the previous rotations. c jm1 = j - 1 if (jm1 .lt. 1) go to 20 do 10 i = 1, jm1 t = c(i)*r(i,j) + s(i)*xj xj = c(i)*xj - s(i)*r(i,j) r(i,j) = t 10 continue 20 continue c c compute the next rotation. c call srotg(r(j,j),xj,c(j),s(j)) 30 continue c c if required, update z and rho. c if (nz .lt. 1) go to 70 do 60 j = 1, nz zeta = y(j) do 40 i = 1, p t = c(i)*z(i,j) + s(i)*zeta zeta = c(i)*zeta - s(i)*z(i,j) z(i,j) = t 40 continue azeta = abs(zeta) if (azeta .eq. 0.0e0 .or. rho(j) .lt. 0.0e0) go to 50 scale = azeta + rho(j) rho(j) = scale*sqrt((azeta/scale)**2+(rho(j)/scale)**2) 50 continue 60 continue 70 continue return end subroutine sgbco(abd,lda,n,ml,mu,ipvt,rcond,z) c*********************************************************************72 integer lda,n,ml,mu,ipvt(1) real abd(lda,1),z(1) real rcond c c sgbco factors a real band matrix by gaussian c elimination and estimates the condition of the matrix. c c if rcond is not needed, sgbfa is slightly faster. c to solve a*x = b , follow sgbco by sgbsl. c to compute inverse(a)*c , follow sgbco by sgbsl. c to compute determinant(a) , follow sgbco by sgbdi. c c on entry c c abd real(lda, n) c contains the matrix in band storage. the columns c of the matrix are stored in the columns of abd and c the diagonals of the matrix are stored in rows c ml+1 through 2*ml+mu+1 of abd . c see the comments below for details. c c lda integer c the leading dimension of the array abd . c lda must be .ge. 2*ml + mu + 1 . c c n integer c the order of the original matrix. c c ml integer c number of diagonals below the main diagonal. c 0 .le. ml .lt. n . c c mu integer c number of diagonals above the main diagonal. c 0 .le. mu .lt. n . c more efficient if ml .le. mu . c c on return c c abd an upper triangular matrix in band storage and c the multipliers which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c c ipvt integer(n) c an integer vector of pivot indices. c c rcond real c an estimate of the reciprocal condition of a . c for the system a*x = b , relative perturbations c in a and b of size epsilon may cause c relative perturbations in x of size epsilon/rcond . c if rcond is so small that the logical expression c 1.0 + rcond .eq. 1.0 c is true, then a may be singular to working c precision. in particular, rcond is zero if c exact singularity is detected or the estimate c underflows. c c z real(n) c a work vector whose contents are usually unimportant. c if a is close to a singular matrix, then z is c an approximate null vector in the sense that c norm(a*z) = rcond*norm(a)*norm(z) . c c band storage c c if a is a band matrix, the following program segment c will set up the input. c c ml = (band width below the diagonal) c mu = (band width above the diagonal) c m = ml + mu + 1 c do 20 j = 1, n c i1 = max0(1, j-mu) c i2 = min0(n, j+ml) c do 10 i = i1, i2 c k = i - j + m c abd(k,j) = a(i,j) c 10 continue c 20 continue c c this uses rows ml+1 through 2*ml+mu+1 of abd . c in addition, the first ml rows in abd are used for c elements generated during the triangularization. c the total number of rows needed in abd is 2*ml+mu+1 . c the ml+mu by ml+mu upper left triangle and the c ml by ml lower right triangle are not referenced. c c example.. if the original matrix is c c 11 12 13 0 0 0 c 21 22 23 24 0 0 c 0 32 33 34 35 0 c 0 0 43 44 45 46 c 0 0 0 54 55 56 c 0 0 0 0 65 66 c c then n = 6, ml = 1, mu = 2, lda .ge. 5 and abd should contain c c * * * + + + , * = not used c * * 13 24 35 46 , + = used for pivoting c * 12 23 34 45 56 c 11 22 33 44 55 66 c 21 32 43 54 65 * c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c linpack sgbfa c blas saxpy,sdot,sscal,sasum c fortran abs,amax1,max0,min0,sign c c internal variables c real sdot,ek,t,wk,wkm real anorm,s,sasum,sm,ynorm integer is,info,j,ju,k,kb,kp1,l,la,lm,lz,m,mm c c c compute 1-norm of a c anorm = 0.0e0 l = ml + 1 is = l + mu do 10 j = 1, n anorm = amax1(anorm,sasum(l,abd(is,j),1)) if (is .gt. ml + 1) is = is - 1 if (j .le. mu) l = l + 1 if (j .ge. n - ml) l = l - 1 10 continue c c factor c call sgbfa(abd,lda,n,ml,mu,ipvt,info) c c rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) . c estimate = norm(z)/norm(y) where a*z = y and trans(a)*y = e . c trans(a) is the transpose of a . the components of e are c chosen to cause maximum local growth in the elements of w where c trans(u)*w = e . the vectors are frequently rescaled to avoid c overflow. c c solve trans(u)*w = e c ek = 1.0e0 do 20 j = 1, n z(j) = 0.0e0 20 continue m = ml + mu + 1 ju = 0 do 100 k = 1, n if (z(k) .ne. 0.0e0) ek = sign(ek,-z(k)) if (abs(ek-z(k)) .le. abs(abd(m,k))) go to 30 s = abs(abd(m,k))/abs(ek-z(k)) call sscal(n,s,z,1) ek = s*ek 30 continue wk = ek - z(k) wkm = -ek - z(k) s = abs(wk) sm = abs(wkm) if (abd(m,k) .eq. 0.0e0) go to 40 wk = wk/abd(m,k) wkm = wkm/abd(m,k) go to 50 40 continue wk = 1.0e0 wkm = 1.0e0 50 continue kp1 = k + 1 ju = min0(max0(ju,mu+ipvt(k)),n) mm = m if (kp1 .gt. ju) go to 90 do 60 j = kp1, ju mm = mm - 1 sm = sm + abs(z(j)+wkm*abd(mm,j)) z(j) = z(j) + wk*abd(mm,j) s = s + abs(z(j)) 60 continue if (s .ge. sm) go to 80 t = wkm - wk wk = wkm mm = m do 70 j = kp1, ju mm = mm - 1 z(j) = z(j) + t*abd(mm,j) 70 continue 80 continue 90 continue z(k) = wk 100 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) c c solve trans(l)*y = w c do 120 kb = 1, n k = n + 1 - kb lm = min0(ml,n-k) if (k .lt. n) z(k) = z(k) + sdot(lm,abd(m+1,k),1,z(k+1),1) if (abs(z(k)) .le. 1.0e0) go to 110 s = 1.0e0/abs(z(k)) call sscal(n,s,z,1) 110 continue l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t 120 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) c ynorm = 1.0e0 c c solve l*v = y c do 140 k = 1, n l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t lm = min0(ml,n-k) if (k .lt. n) call saxpy(lm,t,abd(m+1,k),1,z(k+1),1) if (abs(z(k)) .le. 1.0e0) go to 130 s = 1.0e0/abs(z(k)) call sscal(n,s,z,1) ynorm = s*ynorm 130 continue 140 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) ynorm = s*ynorm c c solve u*z = w c do 160 kb = 1, n k = n + 1 - kb if (abs(z(k)) .le. abs(abd(m,k))) go to 150 s = abs(abd(m,k))/abs(z(k)) call sscal(n,s,z,1) ynorm = s*ynorm 150 continue if (abd(m,k) .ne. 0.0e0) z(k) = z(k)/abd(m,k) if (abd(m,k) .eq. 0.0e0) z(k) = 1.0e0 lm = min0(k,m) - 1 la = m - lm lz = k - lm t = -z(k) call saxpy(lm,t,abd(la,k),1,z(lz),1) 160 continue c make znorm = 1.0 s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) ynorm = s*ynorm c if (anorm .ne. 0.0e0) rcond = ynorm/anorm if (anorm .eq. 0.0e0) rcond = 0.0e0 return end subroutine sgbdi(abd,lda,n,ml,mu,ipvt,det) c*********************************************************************72 integer lda,n,ml,mu,ipvt(1) real abd(lda,1),det(2) c c sgbdi computes the determinant of a band matrix c using the factors computed by sgbco or sgbfa. c if the inverse is needed, use sgbsl n times. c c on entry c c abd real(lda, n) c the output from sgbco or sgbfa. c c lda integer c the leading dimension of the array abd . c c n integer c the order of the original matrix. c c ml integer c number of diagonals below the main diagonal. c c mu integer c number of diagonals above the main diagonal. c c ipvt integer(n) c the pivot vector from sgbco or sgbfa. c c on return c c det real(2) c determinant of original matrix. c determinant = det(1) * 10.0**det(2) c with 1.0 .le. abs(det(1)) .lt. 10.0 c or det(1) = 0.0 . c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c fortran abs c c internal variables c real ten integer i,m c c m = ml + mu + 1 det(1) = 1.0e0 det(2) = 0.0e0 ten = 10.0e0 do 50 i = 1, n if (ipvt(i) .ne. i) det(1) = -det(1) det(1) = abd(m,i)*det(1) c ...exit if (det(1) .eq. 0.0e0) go to 60 10 if (abs(det(1)) .ge. 1.0e0) go to 20 det(1) = ten*det(1) det(2) = det(2) - 1.0e0 go to 10 20 continue 30 if (abs(det(1)) .lt. ten) go to 40 det(1) = det(1)/ten det(2) = det(2) + 1.0e0 go to 30 40 continue 50 continue 60 continue return end subroutine sgbfa(abd,lda,n,ml,mu,ipvt,info) c*********************************************************************72 integer lda,n,ml,mu,ipvt(1),info real abd(lda,1) c c sgbfa factors a real band matrix by elimination. c c sgbfa is usually called by sgbco, but it can be called c directly with a saving in time if rcond is not needed. c c on entry c c abd real(lda, n) c contains the matrix in band storage. the columns c of the matrix are stored in the columns of abd and c the diagonals of the matrix are stored in rows c ml+1 through 2*ml+mu+1 of abd . c see the comments below for details. c c lda integer c the leading dimension of the array abd . c lda must be .ge. 2*ml + mu + 1 . c c n integer c the order of the original matrix. c c ml integer c number of diagonals below the main diagonal. c 0 .le. ml .lt. n . c c mu integer c number of diagonals above the main diagonal. c 0 .le. mu .lt. n . c more efficient if ml .le. mu . c on return c c abd an upper triangular matrix in band storage and c the multipliers which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c c ipvt integer(n) c an integer vector of pivot indices. c c info integer c = 0 normal value. c = k if u(k,k) .eq. 0.0 . this is not an error c condition for this subroutine, but it does c indicate that sgbsl will divide by zero if c called. use rcond in sgbco for a reliable c indication of singularity. c c band storage c c if a is a band matrix, the following program segment c will set up the input. c c ml = (band width below the diagonal) c mu = (band width above the diagonal) c m = ml + mu + 1 c do 20 j = 1, n c i1 = max0(1, j-mu) c i2 = min0(n, j+ml) c do 10 i = i1, i2 c k = i - j + m c abd(k,j) = a(i,j) c 10 continue c 20 continue c c this uses rows ml+1 through 2*ml+mu+1 of abd . c in addition, the first ml rows in abd are used for c elements generated during the triangularization. c the total number of rows needed in abd is 2*ml+mu+1 . c the ml+mu by ml+mu upper left triangle and the c ml by ml lower right triangle are not referenced. c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas saxpy,sscal,isamax c fortran max0,min0 c c internal variables c real t integer i,isamax,i0,j,ju,jz,j0,j1,k,kp1,l,lm,m,mm,nm1 c c m = ml + mu + 1 info = 0 c c zero initial fill-in columns c j0 = mu + 2 j1 = min0(n,m) - 1 if (j1 .lt. j0) go to 30 do 20 jz = j0, j1 i0 = m + 1 - jz do 10 i = i0, ml abd(i,jz) = 0.0e0 10 continue 20 continue 30 continue jz = j1 ju = 0 c c gaussian elimination with partial pivoting c nm1 = n - 1 if (nm1 .lt. 1) go to 130 do 120 k = 1, nm1 kp1 = k + 1 c c zero next fill-in column c jz = jz + 1 if (jz .gt. n) go to 50 if (ml .lt. 1) go to 50 do 40 i = 1, ml abd(i,jz) = 0.0e0 40 continue 50 continue c c find l = pivot index c lm = min0(ml,n-k) l = isamax(lm+1,abd(m,k),1) + m - 1 ipvt(k) = l + k - m c c zero pivot implies this column already triangularized c if (abd(l,k) .eq. 0.0e0) go to 100 c c interchange if necessary c if (l .eq. m) go to 60 t = abd(l,k) abd(l,k) = abd(m,k) abd(m,k) = t 60 continue c c compute multipliers c t = -1.0e0/abd(m,k) call sscal(lm,t,abd(m+1,k),1) c c row elimination with column indexing c ju = min0(max0(ju,mu+ipvt(k)),n) mm = m if (ju .lt. kp1) go to 90 do 80 j = kp1, ju l = l - 1 mm = mm - 1 t = abd(l,j) if (l .eq. mm) go to 70 abd(l,j) = abd(mm,j) abd(mm,j) = t 70 continue call saxpy(lm,t,abd(m+1,k),1,abd(mm+1,j),1) 80 continue 90 continue go to 110 100 continue info = k 110 continue 120 continue 130 continue ipvt(n) = n if (abd(m,n) .eq. 0.0e0) info = n return end subroutine sgbsl(abd,lda,n,ml,mu,ipvt,b,job) c*********************************************************************72 integer lda,n,ml,mu,ipvt(1),job real abd(lda,1),b(1) c c sgbsl solves the real band system c a * x = b or trans(a) * x = b c using the factors computed by sgbco or sgbfa. c c on entry c c abd real(lda, n) c the output from sgbco or sgbfa. c c lda integer c the leading dimension of the array abd . c c n integer c the order of the original matrix. c c ml integer c number of diagonals below the main diagonal. c c mu integer c number of diagonals above the main diagonal. c c ipvt integer(n) c the pivot vector from sgbco or sgbfa. c c b real(n) c the right hand side vector. c c job integer c = 0 to solve a*x = b , c = nonzero to solve trans(a)*x = b , where c trans(a) is the transpose. c c on return c c b the solution vector x . c c error condition c c a division by zero will occur if the input factor contains a c zero on the diagonal. technically this indicates singularity c but it is often caused by improper arguments or improper c setting of lda . it will not occur if the subroutines are c called correctly and if sgbco has set rcond .gt. 0.0 c or sgbfa has set info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call sgbco(abd,lda,n,ml,mu,ipvt,rcond,z) c if (rcond is too small) go to ... c do 10 j = 1, p c call sgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0) c 10 continue c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas saxpy,sdot c fortran min0 c c internal variables c real sdot,t integer k,kb,l,la,lb,lm,m,nm1 c m = mu + ml + 1 nm1 = n - 1 if (job .ne. 0) go to 50 c c job = 0 , solve a * x = b c first solve l*y = b c if (ml .eq. 0) go to 30 if (nm1 .lt. 1) go to 30 do 20 k = 1, nm1 lm = min0(ml,n-k) l = ipvt(k) t = b(l) if (l .eq. k) go to 10 b(l) = b(k) b(k) = t 10 continue call saxpy(lm,t,abd(m+1,k),1,b(k+1),1) 20 continue 30 continue c c now solve u*x = y c do 40 kb = 1, n k = n + 1 - kb b(k) = b(k)/abd(m,k) lm = min0(k,m) - 1 la = m - lm lb = k - lm t = -b(k) call saxpy(lm,t,abd(la,k),1,b(lb),1) 40 continue go to 100 50 continue c c job = nonzero, solve trans(a) * x = b c first solve trans(u)*y = b c do 60 k = 1, n lm = min0(k,m) - 1 la = m - lm lb = k - lm t = sdot(lm,abd(la,k),1,b(lb),1) b(k) = (b(k) - t)/abd(m,k) 60 continue c c now solve trans(l)*x = y c if (ml .eq. 0) go to 90 if (nm1 .lt. 1) go to 90 do 80 kb = 1, nm1 k = n - kb lm = min0(ml,n-k) b(k) = b(k) + sdot(lm,abd(m+1,k),1,b(k+1),1) l = ipvt(k) if (l .eq. k) go to 70 t = b(l) b(l) = b(k) b(k) = t 70 continue 80 continue 90 continue 100 continue return end subroutine sgeco(a,lda,n,ipvt,rcond,z) c*********************************************************************72 integer lda,n,ipvt(1) real a(lda,1),z(1) real rcond c c sgeco factors a real matrix by gaussian elimination c and estimates the condition of the matrix. c c if rcond is not needed, sgefa is slightly faster. c to solve a*x = b , follow sgeco by sgesl. c to compute inverse(a)*c , follow sgeco by sgesl. c to compute determinant(a) , follow sgeco by sgedi. c to compute inverse(a) , follow sgeco by sgedi. c c on entry c c a real(lda, n) c the matrix to be factored. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c on return c c a an upper triangular matrix and the multipliers c which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c c ipvt integer(n) c an integer vector of pivot indices. c c rcond real c an estimate of the reciprocal condition of a . c for the system a*x = b , relative perturbations c in a and b of size epsilon may cause c relative perturbations in x of size epsilon/rcond . c if rcond is so small that the logical expression c 1.0 + rcond .eq. 1.0 c is true, then a may be singular to working c precision. in particular, rcond is zero if c exact singularity is detected or the estimate c underflows. c c z real(n) c a work vector whose contents are usually unimportant. c if a is close to a singular matrix, then z is c an approximate null vector in the sense that c norm(a*z) = rcond*norm(a)*norm(z) . c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c linpack sgefa c blas saxpy,sdot,sscal,sasum c fortran abs,amax1,sign c c internal variables c real sdot,ek,t,wk,wkm real anorm,s,sasum,sm,ynorm integer info,j,k,kb,kp1,l c c c compute 1-norm of a c anorm = 0.0e0 do 10 j = 1, n anorm = amax1(anorm,sasum(n,a(1,j),1)) 10 continue c c factor c call sgefa(a,lda,n,ipvt,info) c c rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) . c estimate = norm(z)/norm(y) where a*z = y and trans(a)*y = e . c trans(a) is the transpose of a . the components of e are c chosen to cause maximum local growth in the elements of w where c trans(u)*w = e . the vectors are frequently rescaled to avoid c overflow. c c solve trans(u)*w = e c ek = 1.0e0 do 20 j = 1, n z(j) = 0.0e0 20 continue do 100 k = 1, n if (z(k) .ne. 0.0e0) ek = sign(ek,-z(k)) if (abs(ek-z(k)) .le. abs(a(k,k))) go to 30 s = abs(a(k,k))/abs(ek-z(k)) call sscal(n,s,z,1) ek = s*ek 30 continue wk = ek - z(k) wkm = -ek - z(k) s = abs(wk) sm = abs(wkm) if (a(k,k) .eq. 0.0e0) go to 40 wk = wk/a(k,k) wkm = wkm/a(k,k) go to 50 40 continue wk = 1.0e0 wkm = 1.0e0 50 continue kp1 = k + 1 if (kp1 .gt. n) go to 90 do 60 j = kp1, n sm = sm + abs(z(j)+wkm*a(k,j)) z(j) = z(j) + wk*a(k,j) s = s + abs(z(j)) 60 continue if (s .ge. sm) go to 80 t = wkm - wk wk = wkm do 70 j = kp1, n z(j) = z(j) + t*a(k,j) 70 continue 80 continue 90 continue z(k) = wk 100 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) c c solve trans(l)*y = w c do 120 kb = 1, n k = n + 1 - kb if (k .lt. n) z(k) = z(k) + sdot(n-k,a(k+1,k),1,z(k+1),1) if (abs(z(k)) .le. 1.0e0) go to 110 s = 1.0e0/abs(z(k)) call sscal(n,s,z,1) 110 continue l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t 120 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) c ynorm = 1.0e0 c c solve l*v = y c do 140 k = 1, n l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t if (k .lt. n) call saxpy(n-k,t,a(k+1,k),1,z(k+1),1) if (abs(z(k)) .le. 1.0e0) go to 130 s = 1.0e0/abs(z(k)) call sscal(n,s,z,1) ynorm = s*ynorm 130 continue 140 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) ynorm = s*ynorm c c solve u*z = v c do 160 kb = 1, n k = n + 1 - kb if (abs(z(k)) .le. abs(a(k,k))) go to 150 s = abs(a(k,k))/abs(z(k)) call sscal(n,s,z,1) ynorm = s*ynorm 150 continue if (a(k,k) .ne. 0.0e0) z(k) = z(k)/a(k,k) if (a(k,k) .eq. 0.0e0) z(k) = 1.0e0 t = -z(k) call saxpy(k-1,t,a(1,k),1,z(1),1) 160 continue c make znorm = 1.0 s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) ynorm = s*ynorm c if (anorm .ne. 0.0e0) rcond = ynorm/anorm if (anorm .eq. 0.0e0) rcond = 0.0e0 return end subroutine sgedi(a,lda,n,ipvt,det,work,job) c*********************************************************************72 integer lda,n,ipvt(1),job real a(lda,1),det(2),work(1) c c sgedi computes the determinant and inverse of a matrix c using the factors computed by sgeco or sgefa. c c on entry c c a real(lda, n) c the output from sgeco or sgefa. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c ipvt integer(n) c the pivot vector from sgeco or sgefa. c c work real(n) c work vector. contents destroyed. c c job integer c = 11 both determinant and inverse. c = 01 inverse only. c = 10 determinant only. c c on return c c a inverse of original matrix if requested. c otherwise unchanged. c c det real(2) c determinant of original matrix if requested. c otherwise not referenced. c determinant = det(1) * 10.0**det(2) c with 1.0 .le. abs(det(1)) .lt. 10.0 c or det(1) .eq. 0.0 . c c error condition c c a division by zero will occur if the input factor contains c a zero on the diagonal and the inverse is requested. c it will not occur if the subroutines are called correctly c and if sgeco has set rcond .gt. 0.0 or sgefa has set c info .eq. 0 . c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas saxpy,sscal,sswap c fortran abs,mod c c internal variables c real t real ten integer i,j,k,kb,kp1,l,nm1 c c c compute determinant c if (job/10 .eq. 0) go to 70 det(1) = 1.0e0 det(2) = 0.0e0 ten = 10.0e0 do 50 i = 1, n if (ipvt(i) .ne. i) det(1) = -det(1) det(1) = a(i,i)*det(1) c ...exit if (det(1) .eq. 0.0e0) go to 60 10 if (abs(det(1)) .ge. 1.0e0) go to 20 det(1) = ten*det(1) det(2) = det(2) - 1.0e0 go to 10 20 continue 30 if (abs(det(1)) .lt. ten) go to 40 det(1) = det(1)/ten det(2) = det(2) + 1.0e0 go to 30 40 continue 50 continue 60 continue 70 continue c c compute inverse(u) c if (mod(job,10) .eq. 0) go to 150 do 100 k = 1, n a(k,k) = 1.0e0/a(k,k) t = -a(k,k) call sscal(k-1,t,a(1,k),1) kp1 = k + 1 if (n .lt. kp1) go to 90 do 80 j = kp1, n t = a(k,j) a(k,j) = 0.0e0 call saxpy(k,t,a(1,k),1,a(1,j),1) 80 continue 90 continue 100 continue c c form inverse(u)*inverse(l) c nm1 = n - 1 if (nm1 .lt. 1) go to 140 do 130 kb = 1, nm1 k = n - kb kp1 = k + 1 do 110 i = kp1, n work(i) = a(i,k) a(i,k) = 0.0e0 110 continue do 120 j = kp1, n t = work(j) call saxpy(n,t,a(1,j),1,a(1,k),1) 120 continue l = ipvt(k) if (l .ne. k) call sswap(n,a(1,k),1,a(1,l),1) 130 continue 140 continue 150 continue return end subroutine sgefa(a,lda,n,ipvt,info) c*********************************************************************72 integer lda,n,ipvt(1),info real a(lda,1) c c sgefa factors a real matrix by gaussian elimination. c c sgefa is usually called by sgeco, but it can be called c directly with a saving in time if rcond is not needed. c (time for sgeco) = (1 + 9/n)*(time for sgefa) . c c on entry c c a real(lda, n) c the matrix to be factored. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c on return c c a an upper triangular matrix and the multipliers c which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c c ipvt integer(n) c an integer vector of pivot indices. c c info integer c = 0 normal value. c = k if u(k,k) .eq. 0.0 . this is not an error c condition for this subroutine, but it does c indicate that sgesl or sgedi will divide by zero c if called. use rcond in sgeco for a reliable c indication of singularity. c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas saxpy,sscal,isamax c c internal variables c real t integer isamax,j,k,kp1,l,nm1 c c c gaussian elimination with partial pivoting c info = 0 nm1 = n - 1 if (nm1 .lt. 1) go to 70 do 60 k = 1, nm1 kp1 = k + 1 c c find l = pivot index c l = isamax(n-k+1,a(k,k),1) + k - 1 ipvt(k) = l c c zero pivot implies this column already triangularized c if (a(l,k) .eq. 0.0e0) go to 40 c c interchange if necessary c if (l .eq. k) go to 10 t = a(l,k) a(l,k) = a(k,k) a(k,k) = t 10 continue c c compute multipliers c t = -1.0e0/a(k,k) call sscal(n-k,t,a(k+1,k),1) c c row elimination with column indexing c do 30 j = kp1, n t = a(l,j) if (l .eq. k) go to 20 a(l,j) = a(k,j) a(k,j) = t 20 continue call saxpy(n-k,t,a(k+1,k),1,a(k+1,j),1) 30 continue go to 50 40 continue info = k 50 continue 60 continue 70 continue ipvt(n) = n if (a(n,n) .eq. 0.0e0) info = n return end subroutine sgesl(a,lda,n,ipvt,b,job) c*********************************************************************72 integer lda,n,ipvt(1),job real a(lda,1),b(1) c c sgesl solves the real system c a * x = b or trans(a) * x = b c using the factors computed by sgeco or sgefa. c c on entry c c a real(lda, n) c the output from sgeco or sgefa. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c ipvt integer(n) c the pivot vector from sgeco or sgefa. c c b real(n) c the right hand side vector. c c job integer c = 0 to solve a*x = b , c = nonzero to solve trans(a)*x = b where c trans(a) is the transpose. c c on return c c b the solution vector x . c c error condition c c a division by zero will occur if the input factor contains a c zero on the diagonal. technically this indicates singularity c but it is often caused by improper arguments or improper c setting of lda . it will not occur if the subroutines are c called correctly and if sgeco has set rcond .gt. 0.0 c or sgefa has set info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call sgeco(a,lda,n,ipvt,rcond,z) c if (rcond is too small) go to ... c do 10 j = 1, p c call sgesl(a,lda,n,ipvt,c(1,j),0) c 10 continue c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas saxpy,sdot c c internal variables c real sdot,t integer k,kb,l,nm1 c nm1 = n - 1 if (job .ne. 0) go to 50 c c job = 0 , solve a * x = b c first solve l*y = b c if (nm1 .lt. 1) go to 30 do 20 k = 1, nm1 l = ipvt(k) t = b(l) if (l .eq. k) go to 10 b(l) = b(k) b(k) = t 10 continue call saxpy(n-k,t,a(k+1,k),1,b(k+1),1) 20 continue 30 continue c c now solve u*x = y c do 40 kb = 1, n k = n + 1 - kb b(k) = b(k)/a(k,k) t = -b(k) call saxpy(k-1,t,a(1,k),1,b(1),1) 40 continue go to 100 50 continue c c job = nonzero, solve trans(a) * x = b c first solve trans(u)*y = b c do 60 k = 1, n t = sdot(k-1,a(1,k),1,b(1),1) b(k) = (b(k) - t)/a(k,k) 60 continue c c now solve trans(l)*x = y c if (nm1 .lt. 1) go to 90 do 80 kb = 1, nm1 k = n - kb b(k) = b(k) + sdot(n-k,a(k+1,k),1,b(k+1),1) l = ipvt(k) if (l .eq. k) go to 70 t = b(l) b(l) = b(k) b(k) = t 70 continue 80 continue 90 continue 100 continue return end subroutine sgtsl(n,c,d,e,b,info) c*********************************************************************72 integer n,info real c(1),d(1),e(1),b(1) c c sgtsl given a general tridiagonal matrix and a right hand c side will find the solution. c c on entry c c n integer c is the order of the tridiagonal matrix. c c c real(n) c is the subdiagonal of the tridiagonal matrix. c c(2) through c(n) should contain the subdiagonal. c on output c is destroyed. c c d real(n) c is the diagonal of the tridiagonal matrix. c on output d is destroyed. c c e real(n) c is the superdiagonal of the tridiagonal matrix. c e(1) through e(n-1) should contain the superdiagonal. c on output e is destroyed. c c b real(n) c is the right hand side vector. c c on return c c b is the solution vector. c c info integer c = 0 normal value. c = k if the k-th element of the diagonal becomes c exactly zero. the subroutine returns when c this is detected. c c linpack. this version dated 08/14/78 . c jack dongarra, argonne national laboratory. c c no externals c fortran abs c c internal variables c integer k,kb,kp1,nm1,nm2 real t c begin block permitting ...exits to 100 c info = 0 c(1) = d(1) nm1 = n - 1 if (nm1 .lt. 1) go to 40 d(1) = e(1) e(1) = 0.0e0 e(n) = 0.0e0 c do 30 k = 1, nm1 kp1 = k + 1 c c find the largest of the two rows c if (abs(c(kp1)) .lt. abs(c(k))) go to 10 c c interchange row c t = c(kp1) c(kp1) = c(k) c(k) = t t = d(kp1) d(kp1) = d(k) d(k) = t t = e(kp1) e(kp1) = e(k) e(k) = t t = b(kp1) b(kp1) = b(k) b(k) = t 10 continue c c zero elements c if (c(k) .ne. 0.0e0) go to 20 info = k c ............exit go to 100 20 continue t = -c(kp1)/c(k) c(kp1) = d(kp1) + t*d(k) d(kp1) = e(kp1) + t*e(k) e(kp1) = 0.0e0 b(kp1) = b(kp1) + t*b(k) 30 continue 40 continue if (c(n) .ne. 0.0e0) go to 50 info = n go to 90 50 continue c c back solve c nm2 = n - 2 b(n) = b(n)/c(n) if (n .eq. 1) go to 80 b(nm1) = (b(nm1) - d(nm1)*b(n))/c(nm1) if (nm2 .lt. 1) go to 70 do 60 kb = 1, nm2 k = nm2 - kb + 1 b(k) = (b(k) - d(k)*b(k+1) - e(k)*b(k+2))/c(k) 60 continue 70 continue 80 continue 90 continue 100 continue return end subroutine spbco(abd,lda,n,m,rcond,z,info) c*********************************************************************72 integer lda,n,m,info real abd(lda,1),z(1) real rcond c c spbco factors a real symmetric positive definite matrix c stored in band form and estimates the condition of the matrix. c c if rcond is not needed, spbfa is slightly faster. c to solve a*x = b , follow spbco by spbsl. c to compute inverse(a)*c , follow spbco by spbsl. c to compute determinant(a) , follow spbco by spbdi. c c on entry c c abd real(lda, n) c the matrix to be factored. the columns of the upper c triangle are stored in the columns of abd and the c diagonals of the upper triangle are stored in the c rows of abd . see the comments below for details. c c lda integer c the leading dimension of the array abd . c lda must be .ge. m + 1 . c c n integer c the order of the matrix a . c c m integer c the number of diagonals above the main diagonal. c 0 .le. m .lt. n . c c on return c c abd an upper triangular matrix r , stored in band c form, so that a = trans(r)*r . c if info .ne. 0 , the factorization is not complete. c c rcond real c an estimate of the reciprocal condition of a . c for the system a*x = b , relative perturbations c in a and b of size epsilon may cause c relative perturbations in x of size epsilon/rcond . c if rcond is so small that the logical expression c 1.0 + rcond .eq. 1.0 c is true, then a may be singular to working c precision. in particular, rcond is zero if c exact singularity is detected or the estimate c underflows. if info .ne. 0 , rcond is unchanged. c c z real(n) c a work vector whose contents are usually unimportant. c if a is singular to working precision, then z is c an approximate null vector in the sense that c norm(a*z) = rcond*norm(a)*norm(z) . c if info .ne. 0 , z is unchanged. c c info integer c = 0 for normal return. c = k signals an error condition. the leading minor c of order k is not positive definite. c c band storage c c if a is a symmetric positive definite band matrix, c the following program segment will set up the input. c c m = (band width above diagonal) c do 20 j = 1, n c i1 = max0(1, j-m) c do 10 i = i1, j c k = i-j+m+1 c abd(k,j) = a(i,j) c 10 continue c 20 continue c c this uses m + 1 rows of a , except for the m by m c upper left triangle, which is ignored. c c example.. if the original matrix is c c 11 12 13 0 0 0 c 12 22 23 24 0 0 c 13 23 33 34 35 0 c 0 24 34 44 45 46 c 0 0 35 45 55 56 c 0 0 0 46 56 66 c c then n = 6 , m = 2 and abd should contain c c * * 13 24 35 46 c * 12 23 34 45 56 c 11 22 33 44 55 66 c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c linpack spbfa c blas saxpy,sdot,sscal,sasum c fortran abs,amax1,max0,min0,real,sign c c internal variables c real sdot,ek,t,wk,wkm real anorm,s,sasum,sm,ynorm integer i,j,j2,k,kb,kp1,l,la,lb,lm,mu c c c find norm of a c do 30 j = 1, n l = min0(j,m+1) mu = max0(m+2-j,1) z(j) = sasum(l,abd(mu,j),1) k = j - l if (m .lt. mu) go to 20 do 10 i = mu, m k = k + 1 z(k) = z(k) + abs(abd(i,j)) 10 continue 20 continue 30 continue anorm = 0.0e0 do 40 j = 1, n anorm = amax1(anorm,z(j)) 40 continue c c factor c call spbfa(abd,lda,n,m,info) if (info .ne. 0) go to 180 c c rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) . c estimate = norm(z)/norm(y) where a*z = y and a*y = e . c the components of e are chosen to cause maximum local c growth in the elements of w where trans(r)*w = e . c the vectors are frequently rescaled to avoid overflow. c c solve trans(r)*w = e c ek = 1.0e0 do 50 j = 1, n z(j) = 0.0e0 50 continue do 110 k = 1, n if (z(k) .ne. 0.0e0) ek = sign(ek,-z(k)) if (abs(ek-z(k)) .le. abd(m+1,k)) go to 60 s = abd(m+1,k)/abs(ek-z(k)) call sscal(n,s,z,1) ek = s*ek 60 continue wk = ek - z(k) wkm = -ek - z(k) s = abs(wk) sm = abs(wkm) wk = wk/abd(m+1,k) wkm = wkm/abd(m+1,k) kp1 = k + 1 j2 = min0(k+m,n) i = m + 1 if (kp1 .gt. j2) go to 100 do 70 j = kp1, j2 i = i - 1 sm = sm + abs(z(j)+wkm*abd(i,j)) z(j) = z(j) + wk*abd(i,j) s = s + abs(z(j)) 70 continue if (s .ge. sm) go to 90 t = wkm - wk wk = wkm i = m + 1 do 80 j = kp1, j2 i = i - 1 z(j) = z(j) + t*abd(i,j) 80 continue 90 continue 100 continue z(k) = wk 110 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) c c solve r*y = w c do 130 kb = 1, n k = n + 1 - kb if (abs(z(k)) .le. abd(m+1,k)) go to 120 s = abd(m+1,k)/abs(z(k)) call sscal(n,s,z,1) 120 continue z(k) = z(k)/abd(m+1,k) lm = min0(k-1,m) la = m + 1 - lm lb = k - lm t = -z(k) call saxpy(lm,t,abd(la,k),1,z(lb),1) 130 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) c ynorm = 1.0e0 c c solve trans(r)*v = y c do 150 k = 1, n lm = min0(k-1,m) la = m + 1 - lm lb = k - lm z(k) = z(k) - sdot(lm,abd(la,k),1,z(lb),1) if (abs(z(k)) .le. abd(m+1,k)) go to 140 s = abd(m+1,k)/abs(z(k)) call sscal(n,s,z,1) ynorm = s*ynorm 140 continue z(k) = z(k)/abd(m+1,k) 150 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) ynorm = s*ynorm c c solve r*z = w c do 170 kb = 1, n k = n + 1 - kb if (abs(z(k)) .le. abd(m+1,k)) go to 160 s = abd(m+1,k)/abs(z(k)) call sscal(n,s,z,1) ynorm = s*ynorm 160 continue z(k) = z(k)/abd(m+1,k) lm = min0(k-1,m) la = m + 1 - lm lb = k - lm t = -z(k) call saxpy(lm,t,abd(la,k),1,z(lb),1) 170 continue c make znorm = 1.0 s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) ynorm = s*ynorm c if (anorm .ne. 0.0e0) rcond = ynorm/anorm if (anorm .eq. 0.0e0) rcond = 0.0e0 180 continue return end subroutine spbdi(abd,lda,n,m,det) c*********************************************************************72 integer lda,n,m real abd(lda,1) real det(2) c c spbdi computes the determinant c of a real symmetric positive definite band matrix c using the factors computed by spbco or spbfa. c if the inverse is needed, use spbsl n times. c c on entry c c abd real(lda, n) c the output from spbco or spbfa. c c lda integer c the leading dimension of the array abd . c c n integer c the order of the matrix a . c c m integer c the number of diagonals above the main diagonal. c c on return c c det real(2) c determinant of original matrix in the form c determinant = det(1) * 10.0**det(2) c with 1.0 .le. det(1) .lt. 10.0 c or det(1) .eq. 0.0 . c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c c internal variables c real s integer i c c compute determinant c det(1) = 1.0e0 det(2) = 0.0e0 s = 10.0e0 do 50 i = 1, n det(1) = abd(m+1,i)**2*det(1) c ...exit if (det(1) .eq. 0.0e0) go to 60 10 if (det(1) .ge. 1.0e0) go to 20 det(1) = s*det(1) det(2) = det(2) - 1.0e0 go to 10 20 continue 30 if (det(1) .lt. s) go to 40 det(1) = det(1)/s det(2) = det(2) + 1.0e0 go to 30 40 continue 50 continue 60 continue return end subroutine spbfa(abd,lda,n,m,info) c*********************************************************************72 integer lda,n,m,info real abd(lda,1) c c spbfa factors a real symmetric positive definite matrix c stored in band form. c c spbfa is usually called by spbco, but it can be called c directly with a saving in time if rcond is not needed. c c on entry c c abd real(lda, n) c the matrix to be factored. the columns of the upper c triangle are stored in the columns of abd and the c diagonals of the upper triangle are stored in the c rows of abd . see the comments below for details. c c lda integer c the leading dimension of the array abd . c lda must be .ge. m + 1 . c c n integer c the order of the matrix a . c c m integer c the number of diagonals above the main diagonal. c 0 .le. m .lt. n . c c on return c c abd an upper triangular matrix r , stored in band c form, so that a = trans(r)*r . c c info integer c = 0 for normal return. c = k if the leading minor of order k is not c positive definite. c c band storage c c if a is a symmetric positive definite band matrix, c the following program segment will set up the input. c c m = (band width above diagonal) c do 20 j = 1, n c i1 = max0(1, j-m) c do 10 i = i1, j c k = i-j+m+1 c abd(k,j) = a(i,j) c 10 continue c 20 continue c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas sdot c fortran max0,sqrt c c internal variables c real sdot,t real s integer ik,j,jk,k,mu c begin block with ...exits to 40 c c do 30 j = 1, n info = j s = 0.0e0 ik = m + 1 jk = max0(j-m,1) mu = max0(m+2-j,1) if (m .lt. mu) go to 20 do 10 k = mu, m t = abd(k,j) - sdot(k-mu,abd(ik,jk),1,abd(mu,j),1) t = t/abd(m+1,jk) abd(k,j) = t s = s + t*t ik = ik - 1 jk = jk + 1 10 continue 20 continue s = abd(m+1,j) - s c ......exit if (s .le. 0.0e0) go to 40 abd(m+1,j) = sqrt(s) 30 continue info = 0 40 continue return end subroutine spbsl(abd,lda,n,m,b) c*********************************************************************72 integer lda,n,m real abd(lda,1),b(1) c c spbsl solves the real symmetric positive definite band c system a*x = b c using the factors computed by spbco or spbfa. c c on entry c c abd real(lda, n) c the output from spbco or spbfa. c c lda integer c the leading dimension of the array abd . c c n integer c the order of the matrix a . c c m integer c the number of diagonals above the main diagonal. c c b real(n) c the right hand side vector. c c on return c c b the solution vector x . c c error condition c c a division by zero will occur if the input factor contains c a zero on the diagonal. technically this indicates c singularity but it is usually caused by improper subroutine c arguments. it will not occur if the subroutines are called c correctly and info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call spbco(abd,lda,n,rcond,z,info) c if (rcond is too small .or. info .ne. 0) go to ... c do 10 j = 1, p c call spbsl(abd,lda,n,c(1,j)) c 10 continue c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas saxpy,sdot c fortran min0 c c internal variables c real sdot,t integer k,kb,la,lb,lm c c solve trans(r)*y = b c do 10 k = 1, n lm = min0(k-1,m) la = m + 1 - lm lb = k - lm t = sdot(lm,abd(la,k),1,b(lb),1) b(k) = (b(k) - t)/abd(m+1,k) 10 continue c c solve r*x = y c do 20 kb = 1, n k = n + 1 - kb lm = min0(k-1,m) la = m + 1 - lm lb = k - lm b(k) = b(k)/abd(m+1,k) t = -b(k) call saxpy(lm,t,abd(la,k),1,b(lb),1) 20 continue return end subroutine spoco(a,lda,n,rcond,z,info) c*********************************************************************72 integer lda,n,info real a(lda,1),z(1) real rcond c c spoco factors a real symmetric positive definite matrix c and estimates the condition of the matrix. c c if rcond is not needed, spofa is slightly faster. c to solve a*x = b , follow spoco by sposl. c to compute inverse(a)*c , follow spoco by sposl. c to compute determinant(a) , follow spoco by spodi. c to compute inverse(a) , follow spoco by spodi. c c on entry c c a real(lda, n) c the symmetric matrix to be factored. only the c diagonal and upper triangle are used. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c on return c c a an upper triangular matrix r so that a = trans(r)*r c where trans(r) is the transpose. c the strict lower triangle is unaltered. c if info .ne. 0 , the factorization is not complete. c c rcond real c an estimate of the reciprocal condition of a . c for the system a*x = b , relative perturbations c in a and b of size epsilon may cause c relative perturbations in x of size epsilon/rcond . c if rcond is so small that the logical expression c 1.0 + rcond .eq. 1.0 c is true, then a may be singular to working c precision. in particular, rcond is zero if c exact singularity is detected or the estimate c underflows. if info .ne. 0 , rcond is unchanged. c c z real(n) c a work vector whose contents are usually unimportant. c if a is close to a singular matrix, then z is c an approximate null vector in the sense that c norm(a*z) = rcond*norm(a)*norm(z) . c if info .ne. 0 , z is unchanged. c c info integer c = 0 for normal return. c = k signals an error condition. the leading minor c of order k is not positive definite. c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c linpack spofa c blas saxpy,sdot,sscal,sasum c fortran abs,amax1,real,sign c c internal variables c real sdot,ek,t,wk,wkm real anorm,s,sasum,sm,ynorm integer i,j,jm1,k,kb,kp1 c c c find norm of a using only upper half c do 30 j = 1, n z(j) = sasum(j,a(1,j),1) jm1 = j - 1 if (jm1 .lt. 1) go to 20 do 10 i = 1, jm1 z(i) = z(i) + abs(a(i,j)) 10 continue 20 continue 30 continue anorm = 0.0e0 do 40 j = 1, n anorm = amax1(anorm,z(j)) 40 continue c c factor c call spofa(a,lda,n,info) if (info .ne. 0) go to 180 c c rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) . c estimate = norm(z)/norm(y) where a*z = y and a*y = e . c the components of e are chosen to cause maximum local c growth in the elements of w where trans(r)*w = e . c the vectors are frequently rescaled to avoid overflow. c c solve trans(r)*w = e c ek = 1.0e0 do 50 j = 1, n z(j) = 0.0e0 50 continue do 110 k = 1, n if (z(k) .ne. 0.0e0) ek = sign(ek,-z(k)) if (abs(ek-z(k)) .le. a(k,k)) go to 60 s = a(k,k)/abs(ek-z(k)) call sscal(n,s,z,1) ek = s*ek 60 continue wk = ek - z(k) wkm = -ek - z(k) s = abs(wk) sm = abs(wkm) wk = wk/a(k,k) wkm = wkm/a(k,k) kp1 = k + 1 if (kp1 .gt. n) go to 100 do 70 j = kp1, n sm = sm + abs(z(j)+wkm*a(k,j)) z(j) = z(j) + wk*a(k,j) s = s + abs(z(j)) 70 continue if (s .ge. sm) go to 90 t = wkm - wk wk = wkm do 80 j = kp1, n z(j) = z(j) + t*a(k,j) 80 continue 90 continue 100 continue z(k) = wk 110 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) c c solve r*y = w c do 130 kb = 1, n k = n + 1 - kb if (abs(z(k)) .le. a(k,k)) go to 120 s = a(k,k)/abs(z(k)) call sscal(n,s,z,1) 120 continue z(k) = z(k)/a(k,k) t = -z(k) call saxpy(k-1,t,a(1,k),1,z(1),1) 130 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) c ynorm = 1.0e0 c c solve trans(r)*v = y c do 150 k = 1, n z(k) = z(k) - sdot(k-1,a(1,k),1,z(1),1) if (abs(z(k)) .le. a(k,k)) go to 140 s = a(k,k)/abs(z(k)) call sscal(n,s,z,1) ynorm = s*ynorm 140 continue z(k) = z(k)/a(k,k) 150 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) ynorm = s*ynorm c c solve r*z = v c do 170 kb = 1, n k = n + 1 - kb if (abs(z(k)) .le. a(k,k)) go to 160 s = a(k,k)/abs(z(k)) call sscal(n,s,z,1) ynorm = s*ynorm 160 continue z(k) = z(k)/a(k,k) t = -z(k) call saxpy(k-1,t,a(1,k),1,z(1),1) 170 continue c make znorm = 1.0 s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) ynorm = s*ynorm c if (anorm .ne. 0.0e0) rcond = ynorm/anorm if (anorm .eq. 0.0e0) rcond = 0.0e0 180 continue return end subroutine spodi(a,lda,n,det,job) c*********************************************************************72 integer lda,n,job real a(lda,1) real det(2) c c spodi computes the determinant and inverse of a certain c real symmetric positive definite matrix (see below) c using the factors computed by spoco, spofa or sqrdc. c c on entry c c a real(lda, n) c the output a from spoco or spofa c or the output x from sqrdc. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c job integer c = 11 both determinant and inverse. c = 01 inverse only. c = 10 determinant only. c c on return c c a if spoco or spofa was used to factor a then c spodi produces the upper half of inverse(a) . c if sqrdc was used to decompose x then c spodi produces the upper half of inverse(trans(x)*x) c where trans(x) is the transpose. c elements of a below the diagonal are unchanged. c if the units digit of job is zero, a is unchanged. c c det real(2) c determinant of a or of trans(x)*x if requested. c otherwise not referenced. c determinant = det(1) * 10.0**det(2) c with 1.0 .le. det(1) .lt. 10.0 c or det(1) .eq. 0.0 . c c error condition c c a division by zero will occur if the input factor contains c a zero on the diagonal and the inverse is requested. c it will not occur if the subroutines are called correctly c and if spoco or spofa has set info .eq. 0 . c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas saxpy,sscal c fortran mod c c internal variables c real t real s integer i,j,jm1,k,kp1 c c compute determinant c if (job/10 .eq. 0) go to 70 det(1) = 1.0e0 det(2) = 0.0e0 s = 10.0e0 do 50 i = 1, n det(1) = a(i,i)**2*det(1) c ...exit if (det(1) .eq. 0.0e0) go to 60 10 if (det(1) .ge. 1.0e0) go to 20 det(1) = s*det(1) det(2) = det(2) - 1.0e0 go to 10 20 continue 30 if (det(1) .lt. s) go to 40 det(1) = det(1)/s det(2) = det(2) + 1.0e0 go to 30 40 continue 50 continue 60 continue 70 continue c c compute inverse(r) c if (mod(job,10) .eq. 0) go to 140 do 100 k = 1, n a(k,k) = 1.0e0/a(k,k) t = -a(k,k) call sscal(k-1,t,a(1,k),1) kp1 = k + 1 if (n .lt. kp1) go to 90 do 80 j = kp1, n t = a(k,j) a(k,j) = 0.0e0 call saxpy(k,t,a(1,k),1,a(1,j),1) 80 continue 90 continue 100 continue c c form inverse(r) * trans(inverse(r)) c do 130 j = 1, n jm1 = j - 1 if (jm1 .lt. 1) go to 120 do 110 k = 1, jm1 t = a(k,j) call saxpy(k,t,a(1,j),1,a(1,k),1) 110 continue 120 continue t = a(j,j) call sscal(j,t,a(1,j),1) 130 continue 140 continue return end subroutine spofa(a,lda,n,info) c*********************************************************************72 integer lda,n,info real a(lda,1) c c spofa factors a real symmetric positive definite matrix. c c spofa is usually called by spoco, but it can be called c directly with a saving in time if rcond is not needed. c (time for spoco) = (1 + 18/n)*(time for spofa) . c c on entry c c a real(lda, n) c the symmetric matrix to be factored. only the c diagonal and upper triangle are used. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c on return c c a an upper triangular matrix r so that a = trans(r)*r c where trans(r) is the transpose. c the strict lower triangle is unaltered. c if info .ne. 0 , the factorization is not complete. c c info integer c = 0 for normal return. c = k signals an error condition. the leading minor c of order k is not positive definite. c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas sdot c fortran sqrt c c internal variables c real sdot,t real s integer j,jm1,k c begin block with ...exits to 40 c c do 30 j = 1, n info = j s = 0.0e0 jm1 = j - 1 if (jm1 .lt. 1) go to 20 do 10 k = 1, jm1 t = a(k,j) - sdot(k-1,a(1,k),1,a(1,j),1) t = t/a(k,k) a(k,j) = t s = s + t*t 10 continue 20 continue s = a(j,j) - s c ......exit if (s .le. 0.0e0) go to 40 a(j,j) = sqrt(s) 30 continue info = 0 40 continue return end subroutine sposl(a,lda,n,b) c*********************************************************************72 integer lda,n real a(lda,1),b(1) c c sposl solves the real symmetric positive definite system c a * x = b c using the factors computed by spoco or spofa. c c on entry c c a real(lda, n) c the output from spoco or spofa. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c b real(n) c the right hand side vector. c c on return c c b the solution vector x . c c error condition c c a division by zero will occur if the input factor contains c a zero on the diagonal. technically this indicates c singularity but it is usually caused by improper subroutine c arguments. it will not occur if the subroutines are called c correctly and info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call spoco(a,lda,n,rcond,z,info) c if (rcond is too small .or. info .ne. 0) go to ... c do 10 j = 1, p c call sposl(a,lda,n,c(1,j)) c 10 continue c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas saxpy,sdot c c internal variables c real sdot,t integer k,kb c c solve trans(r)*y = b c do 10 k = 1, n t = sdot(k-1,a(1,k),1,b(1),1) b(k) = (b(k) - t)/a(k,k) 10 continue c c solve r*x = y c do 20 kb = 1, n k = n + 1 - kb b(k) = b(k)/a(k,k) t = -b(k) call saxpy(k-1,t,a(1,k),1,b(1),1) 20 continue return end subroutine sppco(ap,n,rcond,z,info) c*********************************************************************72 integer n,info real ap(1),z(1) real rcond c c sppco factors a real symmetric positive definite matrix c stored in packed form c and estimates the condition of the matrix. c c if rcond is not needed, sppfa is slightly faster. c to solve a*x = b , follow sppco by sppsl. c to compute inverse(a)*c , follow sppco by sppsl. c to compute determinant(a) , follow sppco by sppdi. c to compute inverse(a) , follow sppco by sppdi. c c on entry c c ap real (n*(n+1)/2) c the packed form of a symmetric matrix a . the c columns of the upper triangle are stored sequentially c in a one-dimensional array of length n*(n+1)/2 . c see comments below for details. c c n integer c the order of the matrix a . c c on return c c ap an upper triangular matrix r , stored in packed c form, so that a = trans(r)*r . c if info .ne. 0 , the factorization is not complete. c c rcond real c an estimate of the reciprocal condition of a . c for the system a*x = b , relative perturbations c in a and b of size epsilon may cause c relative perturbations in x of size epsilon/rcond . c if rcond is so small that the logical expression c 1.0 + rcond .eq. 1.0 c is true, then a may be singular to working c precision. in particular, rcond is zero if c exact singularity is detected or the estimate c underflows. if info .ne. 0 , rcond is unchanged. c c z real(n) c a work vector whose contents are usually unimportant. c if a is singular to working precision, then z is c an approximate null vector in the sense that c norm(a*z) = rcond*norm(a)*norm(z) . c if info .ne. 0 , z is unchanged. c c info integer c = 0 for normal return. c = k signals an error condition. the leading minor c of order k is not positive definite. c c packed storage c c the following program segment will pack the upper c triangle of a symmetric matrix. c c k = 0 c do 20 j = 1, n c do 10 i = 1, j c k = k + 1 c ap(k) = a(i,j) c 10 continue c 20 continue c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c linpack sppfa c blas saxpy,sdot,sscal,sasum c fortran abs,amax1,real,sign c c internal variables c real sdot,ek,t,wk,wkm real anorm,s,sasum,sm,ynorm integer i,ij,j,jm1,j1,k,kb,kj,kk,kp1 c c c find norm of a c j1 = 1 do 30 j = 1, n z(j) = sasum(j,ap(j1),1) ij = j1 j1 = j1 + j jm1 = j - 1 if (jm1 .lt. 1) go to 20 do 10 i = 1, jm1 z(i) = z(i) + abs(ap(ij)) ij = ij + 1 10 continue 20 continue 30 continue anorm = 0.0e0 do 40 j = 1, n anorm = amax1(anorm,z(j)) 40 continue c c factor c call sppfa(ap,n,info) if (info .ne. 0) go to 180 c c rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) . c estimate = norm(z)/norm(y) where a*z = y and a*y = e . c the components of e are chosen to cause maximum local c growth in the elements of w where trans(r)*w = e . c the vectors are frequently rescaled to avoid overflow. c c solve trans(r)*w = e c ek = 1.0e0 do 50 j = 1, n z(j) = 0.0e0 50 continue kk = 0 do 110 k = 1, n kk = kk + k if (z(k) .ne. 0.0e0) ek = sign(ek,-z(k)) if (abs(ek-z(k)) .le. ap(kk)) go to 60 s = ap(kk)/abs(ek-z(k)) call sscal(n,s,z,1) ek = s*ek 60 continue wk = ek - z(k) wkm = -ek - z(k) s = abs(wk) sm = abs(wkm) wk = wk/ap(kk) wkm = wkm/ap(kk) kp1 = k + 1 kj = kk + k if (kp1 .gt. n) go to 100 do 70 j = kp1, n sm = sm + abs(z(j)+wkm*ap(kj)) z(j) = z(j) + wk*ap(kj) s = s + abs(z(j)) kj = kj + j 70 continue if (s .ge. sm) go to 90 t = wkm - wk wk = wkm kj = kk + k do 80 j = kp1, n z(j) = z(j) + t*ap(kj) kj = kj + j 80 continue 90 continue 100 continue z(k) = wk 110 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) c c solve r*y = w c do 130 kb = 1, n k = n + 1 - kb if (abs(z(k)) .le. ap(kk)) go to 120 s = ap(kk)/abs(z(k)) call sscal(n,s,z,1) 120 continue z(k) = z(k)/ap(kk) kk = kk - k t = -z(k) call saxpy(k-1,t,ap(kk+1),1,z(1),1) 130 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) c ynorm = 1.0e0 c c solve trans(r)*v = y c do 150 k = 1, n z(k) = z(k) - sdot(k-1,ap(kk+1),1,z(1),1) kk = kk + k if (abs(z(k)) .le. ap(kk)) go to 140 s = ap(kk)/abs(z(k)) call sscal(n,s,z,1) ynorm = s*ynorm 140 continue z(k) = z(k)/ap(kk) 150 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) ynorm = s*ynorm c c solve r*z = v c do 170 kb = 1, n k = n + 1 - kb if (abs(z(k)) .le. ap(kk)) go to 160 s = ap(kk)/abs(z(k)) call sscal(n,s,z,1) ynorm = s*ynorm 160 continue z(k) = z(k)/ap(kk) kk = kk - k t = -z(k) call saxpy(k-1,t,ap(kk+1),1,z(1),1) 170 continue c make znorm = 1.0 s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) ynorm = s*ynorm c if (anorm .ne. 0.0e0) rcond = ynorm/anorm if (anorm .eq. 0.0e0) rcond = 0.0e0 180 continue return end subroutine sppdi(ap,n,det,job) c*********************************************************************72 integer n,job real ap(1) real det(2) c c sppdi computes the determinant and inverse c of a real symmetric positive definite matrix c using the factors computed by sppco or sppfa . c c on entry c c ap real (n*(n+1)/2) c the output from sppco or sppfa. c c n integer c the order of the matrix a . c c job integer c = 11 both determinant and inverse. c = 01 inverse only. c = 10 determinant only. c c on return c c ap the upper triangular half of the inverse . c the strict lower triangle is unaltered. c c det real(2) c determinant of original matrix if requested. c otherwise not referenced. c determinant = det(1) * 10.0**det(2) c with 1.0 .le. det(1) .lt. 10.0 c or det(1) .eq. 0.0 . c c error condition c c a division by zero will occur if the input factor contains c a zero on the diagonal and the inverse is requested. c it will not occur if the subroutines are called correctly c and if spoco or spofa has set info .eq. 0 . c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas saxpy,sscal c fortran mod c c internal variables c real t real s integer i,ii,j,jj,jm1,j1,k,kj,kk,kp1,k1 c c compute determinant c if (job/10 .eq. 0) go to 70 det(1) = 1.0e0 det(2) = 0.0e0 s = 10.0e0 ii = 0 do 50 i = 1, n ii = ii + i det(1) = ap(ii)**2*det(1) c ...exit if (det(1) .eq. 0.0e0) go to 60 10 if (det(1) .ge. 1.0e0) go to 20 det(1) = s*det(1) det(2) = det(2) - 1.0e0 go to 10 20 continue 30 if (det(1) .lt. s) go to 40 det(1) = det(1)/s det(2) = det(2) + 1.0e0 go to 30 40 continue 50 continue 60 continue 70 continue c c compute inverse(r) c if (mod(job,10) .eq. 0) go to 140 kk = 0 do 100 k = 1, n k1 = kk + 1 kk = kk + k ap(kk) = 1.0e0/ap(kk) t = -ap(kk) call sscal(k-1,t,ap(k1),1) kp1 = k + 1 j1 = kk + 1 kj = kk + k if (n .lt. kp1) go to 90 do 80 j = kp1, n t = ap(kj) ap(kj) = 0.0e0 call saxpy(k,t,ap(k1),1,ap(j1),1) j1 = j1 + j kj = kj + j 80 continue 90 continue 100 continue c c form inverse(r) * trans(inverse(r)) c jj = 0 do 130 j = 1, n j1 = jj + 1 jj = jj + j jm1 = j - 1 k1 = 1 kj = j1 if (jm1 .lt. 1) go to 120 do 110 k = 1, jm1 t = ap(kj) call saxpy(k,t,ap(j1),1,ap(k1),1) k1 = k1 + k kj = kj + 1 110 continue 120 continue t = ap(jj) call sscal(j,t,ap(j1),1) 130 continue 140 continue return end subroutine sppfa(ap,n,info) c*********************************************************************72 integer n,info real ap(1) c c sppfa factors a real symmetric positive definite matrix c stored in packed form. c c sppfa is usually called by sppco, but it can be called c directly with a saving in time if rcond is not needed. c (time for sppco) = (1 + 18/n)*(time for sppfa) . c c on entry c c ap real (n*(n+1)/2) c the packed form of a symmetric matrix a . the c columns of the upper triangle are stored sequentially c in a one-dimensional array of length n*(n+1)/2 . c see comments below for details. c c n integer c the order of the matrix a . c c on return c c ap an upper triangular matrix r , stored in packed c form, so that a = trans(r)*r . c c info integer c = 0 for normal return. c = k if the leading minor of order k is not c positive definite. c c c packed storage c c the following program segment will pack the upper c triangle of a symmetric matrix. c c k = 0 c do 20 j = 1, n c do 10 i = 1, j c k = k + 1 c ap(k) = a(i,j) c 10 continue c 20 continue c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas sdot c fortran sqrt c c internal variables c real sdot,t real s integer j,jj,jm1,k,kj,kk c begin block with ...exits to 40 c c jj = 0 do 30 j = 1, n info = j s = 0.0e0 jm1 = j - 1 kj = jj kk = 0 if (jm1 .lt. 1) go to 20 do 10 k = 1, jm1 kj = kj + 1 t = ap(kj) - sdot(k-1,ap(kk+1),1,ap(jj+1),1) kk = kk + k t = t/ap(kk) ap(kj) = t s = s + t*t 10 continue 20 continue jj = jj + j s = ap(jj) - s c ......exit if (s .le. 0.0e0) go to 40 ap(jj) = sqrt(s) 30 continue info = 0 40 continue return end subroutine sppsl(ap,n,b) c*********************************************************************72 integer n real ap(1),b(1) c c sppsl solves the real symmetric positive definite system c a * x = b c using the factors computed by sppco or sppfa. c c on entry c c ap real (n*(n+1)/2) c the output from sppco or sppfa. c c n integer c the order of the matrix a . c c b real(n) c the right hand side vector. c c on return c c b the solution vector x . c c error condition c c a division by zero will occur if the input factor contains c a zero on the diagonal. technically this indicates c singularity but it is usually caused by improper subroutine c arguments. it will not occur if the subroutines are called c correctly and info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call sppco(ap,n,rcond,z,info) c if (rcond is too small .or. info .ne. 0) go to ... c do 10 j = 1, p c call sppsl(ap,n,c(1,j)) c 10 continue c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas saxpy,sdot c c internal variables c real sdot,t integer k,kb,kk c kk = 0 do 10 k = 1, n t = sdot(k-1,ap(kk+1),1,b(1),1) kk = kk + k b(k) = (b(k) - t)/ap(kk) 10 continue do 20 kb = 1, n k = n + 1 - kb b(k) = b(k)/ap(kk) kk = kk - k t = -b(k) call saxpy(k-1,t,ap(kk+1),1,b(1),1) 20 continue return end subroutine sptsl(n,d,e,b) c*********************************************************************72 integer n real d(n),e(n),b(n) c c sptsl given a positive definite tridiagonal matrix and a right c hand side will find the solution. c c on entry c c n integer c is the order of the tridiagonal matrix. c c d real(n) c is the diagonal of the tridiagonal matrix. c on output d is destroyed. c c e real(n) c is the offdiagonal of the tridiagonal matrix. c e(1) through e(n-1) should contain the c offdiagonal. c c b real(n) c is the right hand side vector. c c on return c c b contains the soultion. c c linpack. this version dated 08/14/78 . c jack dongarra, argonne national laboratory. c c no externals c fortran mod c c internal variables c integer k,kbm1,ke,kf,kp1,nm1,nm1d2 real t1,t2 c c check for 1 x 1 case c if (n .ne. 1) go to 10 b(1) = b(1)/d(1) go to 70 10 continue nm1 = n - 1 nm1d2 = nm1/2 if (n .eq. 2) go to 30 kbm1 = n - 1 c c zero top half of subdiagonal and bottom half of c superdiagonal c do 20 k = 1, nm1d2 t1 = 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*e(kbm1) b(kbm1) = b(kbm1) - t2*b(kbm1+1) kbm1 = kbm1 - 1 20 continue 30 continue kp1 = nm1d2 + 1 c c clean up for possible 2 x 2 block at center c if (mod(n,2) .ne. 0) go to 40 t1 = 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 40 continue c c back solve starting at the center, going towards the top c and bottom c b(kp1) = b(kp1)/d(kp1) if (n .eq. 2) go to 60 k = kp1 - 1 ke = kp1 + nm1d2 - 1 do 50 kf = kp1, ke b(k) = (b(k) - e(k)*b(k+1))/d(k) b(kf+1) = (b(kf+1) - e(kf)*b(kf))/d(kf+1) k = k - 1 50 continue 60 continue if (mod(n,2) .eq. 0) b(1) = (b(1) - e(1)*b(2))/d(1) 70 continue return end subroutine sqrdc(x,ldx,n,p,qraux,jpvt,work,job) c*********************************************************************72 integer ldx,n,p,job integer jpvt(1) real x(ldx,1),qraux(1),work(1) c c sqrdc uses householder transformations to compute the qr c factorization of an n by p matrix x. column pivoting c based on the 2-norms of the reduced columns may be c performed at the users option. c c on entry c c x real(ldx,p), where ldx .ge. n. c x contains the matrix whose decomposition is to be c computed. c c ldx integer. c ldx is the leading dimension of the array x. c c n integer. c n is the number of rows of the matrix x. c c p integer. c p is the number of columns of the matrix x. c c jpvt integer(p). c jpvt contains integers that control the selection c of the pivot columns. the k-th column x(k) of x c is placed in one of three classes according to the c value of jpvt(k). c c if jpvt(k) .gt. 0, then x(k) is an initial c column. c c if jpvt(k) .eq. 0, then x(k) is a free column. c c if jpvt(k) .lt. 0, then x(k) is a final column. c c before the decomposition is computed, initial columns c are moved to the beginning of the array x and final c columns to the end. both initial and final columns c are frozen in place during the computation and only c free columns are moved. at the k-th stage of the c reduction, if x(k) is occupied by a free column c it is interchanged with the free column of largest c reduced norm. jpvt is not referenced if c job .eq. 0. c c work real(p). c work is a work array. work is not referenced if c job .eq. 0. c c job integer. c job is an integer that initiates column pivoting. c if job .eq. 0, no pivoting is done. c if job .ne. 0, pivoting is done. c c on return c c x x contains in its upper triangle the upper c triangular matrix r of the qr factorization. c below its diagonal x contains information from c which the orthogonal part of the decomposition c can be recovered. note that if pivoting has c been requested, the decomposition is not that c of the original matrix x but that of x c with its columns permuted as described by jpvt. c c qraux real(p). c qraux contains further information required to recover c the orthogonal part of the decomposition. c c jpvt jpvt(k) contains the index of the column of the c original matrix that has been interchanged into c the k-th column, if pivoting was requested. c c linpack. this version dated 08/14/78 . c g.w. stewart, university of maryland, argonne national lab. c c sqrdc uses the following functions and subprograms. c c blas saxpy,sdot,sscal,sswap,snrm2 c fortran abs,amax1,min0,sqrt c c internal variables c integer j,jp,l,lp1,lup,maxj,pl,pu real maxnrm,snrm2,tt real sdot,nrmxl,t logical negj,swapj c c pl = 1 pu = 0 if (job .eq. 0) go to 60 c c pivoting has been requested. rearrange the columns c according to jpvt. c do 20 j = 1, p swapj = jpvt(j) .gt. 0 negj = jpvt(j) .lt. 0 jpvt(j) = j if (negj) jpvt(j) = -j if (.not.swapj) go to 10 if (j .ne. pl) call sswap(n,x(1,pl),1,x(1,j),1) jpvt(j) = jpvt(pl) jpvt(pl) = j pl = pl + 1 10 continue 20 continue pu = p do 50 jj = 1, p j = p - jj + 1 if (jpvt(j) .ge. 0) go to 40 jpvt(j) = -jpvt(j) if (j .eq. pu) go to 30 call sswap(n,x(1,pu),1,x(1,j),1) jp = jpvt(pu) jpvt(pu) = jpvt(j) jpvt(j) = jp 30 continue pu = pu - 1 40 continue 50 continue 60 continue c c compute the norms of the free columns. c if (pu .lt. pl) go to 80 do 70 j = pl, pu qraux(j) = snrm2(n,x(1,j),1) work(j) = qraux(j) 70 continue 80 continue c c perform the householder reduction of x. c lup = min0(n,p) do 200 l = 1, lup if (l .lt. pl .or. l .ge. pu) go to 120 c c locate the column of largest norm and bring it c into the pivot position. c maxnrm = 0.0e0 maxj = l do 100 j = l, pu if (qraux(j) .le. maxnrm) go to 90 maxnrm = qraux(j) maxj = j 90 continue 100 continue if (maxj .eq. l) go to 110 call sswap(n,x(1,l),1,x(1,maxj),1) qraux(maxj) = qraux(l) work(maxj) = work(l) jp = jpvt(maxj) jpvt(maxj) = jpvt(l) jpvt(l) = jp 110 continue 120 continue qraux(l) = 0.0e0 if (l .eq. n) go to 190 c c compute the householder transformation for column l. c nrmxl = snrm2(n-l+1,x(l,l),1) if (nrmxl .eq. 0.0e0) go to 180 if (x(l,l) .ne. 0.0e0) nrmxl = sign(nrmxl,x(l,l)) call sscal(n-l+1,1.0e0/nrmxl,x(l,l),1) x(l,l) = 1.0e0 + x(l,l) c c apply the transformation to the remaining columns, c updating the norms. c lp1 = l + 1 if (p .lt. lp1) go to 170 do 160 j = lp1, p t = -sdot(n-l+1,x(l,l),1,x(l,j),1)/x(l,l) call saxpy(n-l+1,t,x(l,l),1,x(l,j),1) if (j .lt. pl .or. j .gt. pu) go to 150 if (qraux(j) .eq. 0.0e0) go to 150 tt = 1.0e0 - (abs(x(l,j))/qraux(j))**2 tt = amax1(tt,0.0e0) t = tt tt = 1.0e0 + 0.05e0*tt*(qraux(j)/work(j))**2 if (tt .eq. 1.0e0) go to 130 qraux(j) = qraux(j)*sqrt(t) go to 140 130 continue qraux(j) = snrm2(n-l,x(l+1,j),1) work(j) = qraux(j) 140 continue 150 continue 160 continue 170 continue c c save the transformation. c qraux(l) = x(l,l) x(l,l) = -nrmxl 180 continue 190 continue 200 continue return end subroutine sqrsl(x,ldx,n,k,qraux,y,qy,qty,b,rsd,xb,job,info) c*********************************************************************72 integer ldx,n,k,job,info real x(ldx,1),qraux(1),y(1),qy(1),qty(1),b(1),rsd(1),xb(1) c c sqrsl applies the output of sqrdc to compute coordinate c transformations, projections, and least squares solutions. c for k .le. min(n,p), let xk be the matrix c c xk = (x(jpvt(1)),x(jpvt(2)), ... ,x(jpvt(k))) c c formed from columnns jpvt(1), ... ,jpvt(k) of the original c n x p matrix x that was input to sqrdc (if no pivoting was c done, xk consists of the first k columns of x in their c original order). sqrdc produces a factored orthogonal matrix q c and an upper triangular matrix r such that c c xk = q * (r) c (0) c c this information is contained in coded form in the arrays c x and qraux. c c on entry c c x real(ldx,p). c x contains the output of sqrdc. c c ldx integer. c ldx is the leading dimension of the array x. c c n integer. c n is the number of rows of the matrix xk. it must c have the same value as n in sqrdc. c c k integer. c k is the number of columns of the matrix xk. k c must nnot be greater than min(n,p), where p is the c same as in the calling sequence to sqrdc. c c qraux real(p). c qraux contains the auxiliary output from sqrdc. c c y real(n) c y contains an n-vector that is to be manipulated c by sqrsl. c c job integer. c job specifies what is to be computed. job has c the decimal expansion abcde, with the following c meaning. c c if a.ne.0, compute qy. c if b,c,d, or e .ne. 0, compute qty. c if c.ne.0, compute b. c if d.ne.0, compute rsd. c if e.ne.0, compute xb. c c note that a request to compute b, rsd, or xb c automatically triggers the computation of qty, for c which an array must be provided in the calling c sequence. c c on return c c qy real(n). c qy conntains q*y, if its computation has been c requested. c c qty real(n). c qty contains trans(q)*y, if its computation has c been requested. here trans(q) is the c transpose of the matrix q. c c b real(k) c b contains the solution of the least squares problem c c minimize norm2(y - xk*b), c c if its computation has been requested. (note that c if pivoting was requested in sqrdc, the j-th c component of b will be associated with column jpvt(j) c of the original matrix x that was input into sqrdc.) c c rsd real(n). c rsd contains the least squares residual y - xk*b, c if its computation has been requested. rsd is c also the orthogonal projection of y onto the c orthogonal complement of the column space of xk. c c xb real(n). c xb contains the least squares approximation xk*b, c if its computation has been requested. xb is also c the orthogonal projection of y onto the column space c of x. c c info integer. c info is zero unless the computation of b has c been requested and r is exactly singular. in c this case, info is the index of the first zero c diagonal element of r and b is left unaltered. c c the parameters qy, qty, b, rsd, and xb are not referenced c if their computation is not requested and in this case c can be replaced by dummy variables in the calling program. c to save storage, the user may in some cases use the same c array for different parameters in the calling sequence. a c frequently occuring example is when one wishes to compute c any of b, rsd, or xb and does not need y or qty. in this c case one may identify y, qty, and one of b, rsd, or xb, while c providing separate arrays for anything else that is to be c computed. thus the calling sequence c c call sqrsl(x,ldx,n,k,qraux,y,dum,y,b,y,dum,110,info) c c will result in the computation of b and rsd, with rsd c overwriting y. more generally, each item in the following c list contains groups of permissible identifications for c a single callinng sequence. c c 1. (y,qty,b) (rsd) (xb) (qy) c c 2. (y,qty,rsd) (b) (xb) (qy) c c 3. (y,qty,xb) (b) (rsd) (qy) c c 4. (y,qy) (qty,b) (rsd) (xb) c c 5. (y,qy) (qty,rsd) (b) (xb) c c 6. (y,qy) (qty,xb) (b) (rsd) c c in any group the value returned in the array allocated to c the group corresponds to the last member of the group. c c linpack. this version dated 08/14/78 . c g.w. stewart, university of maryland, argonne national lab. c c sqrsl uses the following functions and subprograms. c c blas saxpy,scopy,sdot c fortran abs,min0,mod c c internal variables c integer i,j,jj,ju,kp1 real sdot,t,temp logical cb,cqy,cqty,cr,cxb c c c set info flag. c info = 0 c c determine what is to be computed. c cqy = job/10000 .ne. 0 cqty = mod(job,10000) .ne. 0 cb = mod(job,1000)/100 .ne. 0 cr = mod(job,100)/10 .ne. 0 cxb = mod(job,10) .ne. 0 ju = min0(k,n-1) c c special action when n=1. c if (ju .ne. 0) go to 40 if (cqy) qy(1) = y(1) if (cqty) qty(1) = y(1) if (cxb) xb(1) = y(1) if (.not.cb) go to 30 if (x(1,1) .ne. 0.0e0) go to 10 info = 1 go to 20 10 continue b(1) = y(1)/x(1,1) 20 continue 30 continue if (cr) rsd(1) = 0.0e0 go to 250 40 continue c c set up to compute qy or qty. c if (cqy) call scopy(n,y,1,qy,1) if (cqty) call scopy(n,y,1,qty,1) if (.not.cqy) go to 70 c c compute qy. c do 60 jj = 1, ju j = ju - jj + 1 if (qraux(j) .eq. 0.0e0) go to 50 temp = x(j,j) x(j,j) = qraux(j) t = -sdot(n-j+1,x(j,j),1,qy(j),1)/x(j,j) call saxpy(n-j+1,t,x(j,j),1,qy(j),1) x(j,j) = temp 50 continue 60 continue 70 continue if (.not.cqty) go to 100 c c compute trans(q)*y. c do 90 j = 1, ju if (qraux(j) .eq. 0.0e0) go to 80 temp = x(j,j) x(j,j) = qraux(j) t = -sdot(n-j+1,x(j,j),1,qty(j),1)/x(j,j) call saxpy(n-j+1,t,x(j,j),1,qty(j),1) x(j,j) = temp 80 continue 90 continue 100 continue c c set up to compute b, rsd, or xb. c if (cb) call scopy(k,qty,1,b,1) kp1 = k + 1 if (cxb) call scopy(k,qty,1,xb,1) if (cr .and. k .lt. n) call scopy(n-k,qty(kp1),1,rsd(kp1),1) if (.not.cxb .or. kp1 .gt. n) go to 120 do 110 i = kp1, n xb(i) = 0.0e0 110 continue 120 continue if (.not.cr) go to 140 do 130 i = 1, k rsd(i) = 0.0e0 130 continue 140 continue if (.not.cb) go to 190 c c compute b. c do 170 jj = 1, k j = k - jj + 1 if (x(j,j) .ne. 0.0e0) go to 150 info = j c ......exit go to 180 150 continue b(j) = b(j)/x(j,j) if (j .eq. 1) go to 160 t = -b(j) call saxpy(j-1,t,x(1,j),1,b,1) 160 continue 170 continue 180 continue 190 continue if (.not.cr .and. .not.cxb) go to 240 c c compute rsd or xb as required. c do 230 jj = 1, ju j = ju - jj + 1 if (qraux(j) .eq. 0.0e0) go to 220 temp = x(j,j) x(j,j) = qraux(j) if (.not.cr) go to 200 t = -sdot(n-j+1,x(j,j),1,rsd(j),1)/x(j,j) call saxpy(n-j+1,t,x(j,j),1,rsd(j),1) 200 continue if (.not.cxb) go to 210 t = -sdot(n-j+1,x(j,j),1,xb(j),1)/x(j,j) call saxpy(n-j+1,t,x(j,j),1,xb(j),1) 210 continue x(j,j) = temp 220 continue 230 continue 240 continue 250 continue return end subroutine ssico(a,lda,n,kpvt,rcond,z) c*********************************************************************72 integer lda,n,kpvt(1) real a(lda,1),z(1) real rcond c c ssico factors a real symmetric matrix by elimination with c symmetric pivoting and estimates the condition of the matrix. c c if rcond is not needed, ssifa is slightly faster. c to solve a*x = b , follow ssico by ssisl. c to compute inverse(a)*c , follow ssico by ssisl. c to compute inverse(a) , follow ssico by ssidi. c to compute determinant(a) , follow ssico by ssidi. c to compute inertia(a), follow ssico by ssidi. c c on entry c c a real(lda, n) c the symmetric matrix to be factored. c only the diagonal and upper triangle are used. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c output c c a a block diagonal matrix and the multipliers which c were used to obtain it. c the factorization can be written a = u*d*trans(u) c where u is a product of permutation and unit c upper triangular matrices , trans(u) is the c transpose of u , and d is block diagonal c with 1 by 1 and 2 by 2 blocks. c c kpvt integer(n) c an integer vector of pivot indices. c c rcond real c an estimate of the reciprocal condition of a . c for the system a*x = b , relative perturbations c in a and b of size epsilon may cause c relative perturbations in x of size epsilon/rcond . c if rcond is so small that the logical expression c 1.0 + rcond .eq. 1.0 c is true, then a may be singular to working c precision. in particular, rcond is zero if c exact singularity is detected or the estimate c underflows. c c z real(n) c a work vector whose contents are usually unimportant. c if a is close to a singular matrix, then z is c an approximate null vector in the sense that c norm(a*z) = rcond*norm(a)*norm(z) . c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c linpack ssifa c blas saxpy,sdot,sscal,sasum c fortran abs,amax1,iabs,sign c c internal variables c real ak,akm1,bk,bkm1,sdot,denom,ek,t real anorm,s,sasum,ynorm integer i,info,j,jm1,k,kp,kps,ks c c c find norm of a using only upper half c do 30 j = 1, n z(j) = sasum(j,a(1,j),1) jm1 = j - 1 if (jm1 .lt. 1) go to 20 do 10 i = 1, jm1 z(i) = z(i) + abs(a(i,j)) 10 continue 20 continue 30 continue anorm = 0.0e0 do 40 j = 1, n anorm = amax1(anorm,z(j)) 40 continue c c factor c call ssifa(a,lda,n,kpvt,info) c c rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) . c estimate = norm(z)/norm(y) where a*z = y and a*y = e . c the components of e are chosen to cause maximum local c growth in the elements of w where u*d*w = e . c the vectors are frequently rescaled to avoid overflow. c c solve u*d*w = e c ek = 1.0e0 do 50 j = 1, n z(j) = 0.0e0 50 continue k = n 60 if (k .eq. 0) go to 120 ks = 1 if (kpvt(k) .lt. 0) ks = 2 kp = iabs(kpvt(k)) kps = k + 1 - ks if (kp .eq. kps) go to 70 t = z(kps) z(kps) = z(kp) z(kp) = t 70 continue if (z(k) .ne. 0.0e0) ek = sign(ek,z(k)) z(k) = z(k) + ek call saxpy(k-ks,z(k),a(1,k),1,z(1),1) if (ks .eq. 1) go to 80 if (z(k-1) .ne. 0.0e0) ek = sign(ek,z(k-1)) z(k-1) = z(k-1) + ek call saxpy(k-ks,z(k-1),a(1,k-1),1,z(1),1) 80 continue if (ks .eq. 2) go to 100 if (abs(z(k)) .le. abs(a(k,k))) go to 90 s = abs(a(k,k))/abs(z(k)) call sscal(n,s,z,1) ek = s*ek 90 continue if (a(k,k) .ne. 0.0e0) z(k) = z(k)/a(k,k) if (a(k,k) .eq. 0.0e0) z(k) = 1.0e0 go to 110 100 continue 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 110 continue k = k - ks go to 60 120 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) c c solve trans(u)*y = w c k = 1 130 if (k .gt. n) go to 160 ks = 1 if (kpvt(k) .lt. 0) ks = 2 if (k .eq. 1) go to 150 z(k) = z(k) + sdot(k-1,a(1,k),1,z(1),1) if (ks .eq. 2) * z(k+1) = z(k+1) + sdot(k-1,a(1,k+1),1,z(1),1) kp = iabs(kpvt(k)) if (kp .eq. k) go to 140 t = z(k) z(k) = z(kp) z(kp) = t 140 continue 150 continue k = k + ks go to 130 160 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) c ynorm = 1.0e0 c c solve u*d*v = y c k = n 170 if (k .eq. 0) go to 230 ks = 1 if (kpvt(k) .lt. 0) ks = 2 if (k .eq. ks) go to 190 kp = iabs(kpvt(k)) kps = k + 1 - ks if (kp .eq. kps) go to 180 t = z(kps) z(kps) = z(kp) z(kp) = t 180 continue call saxpy(k-ks,z(k),a(1,k),1,z(1),1) if (ks .eq. 2) call saxpy(k-ks,z(k-1),a(1,k-1),1,z(1),1) 190 continue if (ks .eq. 2) go to 210 if (abs(z(k)) .le. abs(a(k,k))) go to 200 s = abs(a(k,k))/abs(z(k)) call sscal(n,s,z,1) ynorm = s*ynorm 200 continue if (a(k,k) .ne. 0.0e0) z(k) = z(k)/a(k,k) if (a(k,k) .eq. 0.0e0) z(k) = 1.0e0 go to 220 210 continue 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 220 continue k = k - ks go to 170 230 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) ynorm = s*ynorm c c solve trans(u)*z = v c k = 1 240 if (k .gt. n) go to 270 ks = 1 if (kpvt(k) .lt. 0) ks = 2 if (k .eq. 1) go to 260 z(k) = z(k) + sdot(k-1,a(1,k),1,z(1),1) if (ks .eq. 2) * z(k+1) = z(k+1) + sdot(k-1,a(1,k+1),1,z(1),1) kp = iabs(kpvt(k)) if (kp .eq. k) go to 250 t = z(k) z(k) = z(kp) z(kp) = t 250 continue 260 continue k = k + ks go to 240 270 continue c c make znorm = 1.0 c s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) ynorm = s*ynorm if (anorm .ne. 0.0e0) rcond = ynorm/anorm if (anorm .eq. 0.0e0) rcond = 0.0e0 return end subroutine ssidi(a,lda,n,kpvt,det,inert,work,job) c*********************************************************************72 integer lda,n,job real a(lda,1),work(1) real det(2) integer kpvt(1),inert(3) c c ssidi computes the determinant, inertia and inverse c of a real symmetric matrix using the factors from ssifa. c c on entry c c a real(lda,n) c the output from ssifa. c c lda integer c the leading dimension of the array a. c c n integer c the order of the matrix a. c c kpvt integer(n) c the pivot vector from ssifa. c c work real(n) c work vector. contents destroyed. c c job integer c job has the decimal expansion abc where c if c .ne. 0, the inverse is computed, c if b .ne. 0, the determinant is computed, c if a .ne. 0, the inertia is computed. c c for example, job = 111 gives all three. c c on return c c variables not requested by job are not used. c c a contains the upper triangle of the inverse of c the original matrix. the strict lower triangle c is never referenced. c c det real(2) c determinant of original matrix. c determinant = det(1) * 10.0**det(2) c with 1.0 .le. abs(det(1)) .lt. 10.0 c or det(1) = 0.0. c c inert integer(3) c the inertia of the original matrix. c inert(1) = number of positive eigenvalues. c inert(2) = number of negative eigenvalues. c inert(3) = number of zero eigenvalues. c c error condition c c a division by zero may occur if the inverse is requested c and ssico has set rcond .eq. 0.0 c or ssifa has set info .ne. 0 . c c linpack. this version dated 08/14/78 . c james bunch, univ. calif. san diego, argonne nat. lab c c subroutines and functions c c blas saxpy,scopy,sdot,sswap c fortran abs,iabs,mod c c internal variables. c real akkp1,sdot,temp real ten,d,t,ak,akp1 integer j,jb,k,km1,ks,kstep logical noinv,nodet,noert c noinv = mod(job,10) .eq. 0 nodet = mod(job,100)/10 .eq. 0 noert = mod(job,1000)/100 .eq. 0 c if (nodet .and. noert) go to 140 if (noert) go to 10 inert(1) = 0 inert(2) = 0 inert(3) = 0 10 continue if (nodet) go to 20 det(1) = 1.0e0 det(2) = 0.0e0 ten = 10.0e0 20 continue t = 0.0e0 do 130 k = 1, n d = a(k,k) c c check if 1 by 1 c if (kpvt(k) .gt. 0) go to 50 c c 2 by 2 block c use det (d s) = (d/t * c - t) * t , t = abs(s) c (s c) c to avoid underflow/overflow troubles. c take two passes through scaling. use t for flag. c if (t .ne. 0.0e0) go to 30 t = abs(a(k,k+1)) d = (d/t)*a(k+1,k+1) - t go to 40 30 continue d = t t = 0.0e0 40 continue 50 continue c if (noert) go to 60 if (d .gt. 0.0e0) inert(1) = inert(1) + 1 if (d .lt. 0.0e0) inert(2) = inert(2) + 1 if (d .eq. 0.0e0) inert(3) = inert(3) + 1 60 continue c if (nodet) go to 120 det(1) = d*det(1) if (det(1) .eq. 0.0e0) go to 110 70 if (abs(det(1)) .ge. 1.0e0) go to 80 det(1) = ten*det(1) det(2) = det(2) - 1.0e0 go to 70 80 continue 90 if (abs(det(1)) .lt. ten) go to 100 det(1) = det(1)/ten det(2) = det(2) + 1.0e0 go to 90 100 continue 110 continue 120 continue 130 continue 140 continue c c compute inverse(a) c if (noinv) go to 270 k = 1 150 if (k .gt. n) go to 260 km1 = k - 1 if (kpvt(k) .lt. 0) go to 180 c c 1 by 1 c a(k,k) = 1.0e0/a(k,k) if (km1 .lt. 1) go to 170 call scopy(km1,a(1,k),1,work,1) do 160 j = 1, km1 a(j,k) = sdot(j,a(1,j),1,work,1) call saxpy(j-1,work(j),a(1,j),1,a(1,k),1) 160 continue a(k,k) = a(k,k) + sdot(km1,work,1,a(1,k),1) 170 continue kstep = 1 go to 220 180 continue c c 2 by 2 c t = abs(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) a(k,k) = akp1/d a(k+1,k+1) = ak/d a(k,k+1) = -akkp1/d if (km1 .lt. 1) go to 210 call scopy(km1,a(1,k+1),1,work,1) do 190 j = 1, km1 a(j,k+1) = sdot(j,a(1,j),1,work,1) call saxpy(j-1,work(j),a(1,j),1,a(1,k+1),1) 190 continue a(k+1,k+1) = a(k+1,k+1) + sdot(km1,work,1,a(1,k+1),1) a(k,k+1) = a(k,k+1) + sdot(km1,a(1,k),1,a(1,k+1),1) call scopy(km1,a(1,k),1,work,1) do 200 j = 1, km1 a(j,k) = sdot(j,a(1,j),1,work,1) call saxpy(j-1,work(j),a(1,j),1,a(1,k),1) 200 continue a(k,k) = a(k,k) + sdot(km1,work,1,a(1,k),1) 210 continue kstep = 2 220 continue c c swap c ks = iabs(kpvt(k)) if (ks .eq. k) go to 250 call sswap(ks,a(1,ks),1,a(1,k),1) do 230 jb = ks, k j = k + ks - jb temp = a(j,k) a(j,k) = a(ks,j) a(ks,j) = temp 230 continue if (kstep .eq. 1) go to 240 temp = a(ks,k+1) a(ks,k+1) = a(k,k+1) a(k,k+1) = temp 240 continue 250 continue k = k + kstep go to 150 260 continue 270 continue return end subroutine ssifa(a,lda,n,kpvt,info) c*********************************************************************72 integer lda,n,kpvt(1),info real a(lda,1) c c ssifa factors a real symmetric matrix by elimination c with symmetric pivoting. c c to solve a*x = b , follow ssifa by ssisl. c to compute inverse(a)*c , follow ssifa by ssisl. c to compute determinant(a) , follow ssifa by ssidi. c to compute inertia(a) , follow ssifa by ssidi. c to compute inverse(a) , follow ssifa by ssidi. c c on entry c c a real(lda,n) c the symmetric matrix to be factored. c only the diagonal and upper triangle are used. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c on return c c a a block diagonal matrix and the multipliers which c were used to obtain it. c the factorization can be written a = u*d*trans(u) c where u is a product of permutation and unit c upper triangular matrices , trans(u) is the c transpose of u , and d is block diagonal c with 1 by 1 and 2 by 2 blocks. c c kpvt integer(n) c an integer vector of pivot indices. c c info integer c = 0 normal value. c = k if the k-th pivot block is singular. this is c not an error condition for this subroutine, c but it does indicate that ssisl or ssidi may c divide by zero if called. c c linpack. this version dated 08/14/78 . c james bunch, univ. calif. san diego, argonne nat. lab. c c subroutines and functions c c blas saxpy,sswap,isamax c fortran abs,amax1,sqrt c c internal variables c real ak,akm1,bk,bkm1,denom,mulk,mulkm1,t real absakk,alpha,colmax,rowmax integer imax,imaxp1,j,jj,jmax,k,km1,km2,kstep,isamax logical swap c c c initialize c c alpha is used in choosing pivot block size. alpha = (1.0e0 + sqrt(17.0e0))/8.0e0 c info = 0 c c main loop on k, which goes from n to 1. c k = n 10 continue c c leave the loop if k=0 or k=1. c c ...exit if (k .eq. 0) go to 200 if (k .gt. 1) go to 20 kpvt(1) = 1 if (a(1,1) .eq. 0.0e0) info = 1 c ......exit go to 200 20 continue c c this section of code determines the kind of c elimination to be performed. when it is completed, c kstep will be set to the size of the pivot block, and c swap will be set to .true. if an interchange is c required. c km1 = k - 1 absakk = abs(a(k,k)) c c determine the largest off-diagonal element in c column k. c imax = isamax(k-1,a(1,k),1) colmax = abs(a(imax,k)) if (absakk .lt. alpha*colmax) go to 30 kstep = 1 swap = .false. go to 90 30 continue c c determine the largest off-diagonal element in c row imax. c rowmax = 0.0e0 imaxp1 = imax + 1 do 40 j = imaxp1, k rowmax = amax1(rowmax,abs(a(imax,j))) 40 continue if (imax .eq. 1) go to 50 jmax = isamax(imax-1,a(1,imax),1) rowmax = amax1(rowmax,abs(a(jmax,imax))) 50 continue if (abs(a(imax,imax)) .lt. alpha*rowmax) go to 60 kstep = 1 swap = .true. go to 80 60 continue if (absakk .lt. alpha*colmax*(colmax/rowmax)) go to 70 kstep = 1 swap = .false. go to 80 70 continue kstep = 2 swap = imax .ne. km1 80 continue 90 continue if (amax1(absakk,colmax) .ne. 0.0e0) go to 100 c c column k is zero. set info and iterate the loop. c kpvt(k) = k info = k go to 190 100 continue if (kstep .eq. 2) go to 140 c c 1 x 1 pivot block. c if (.not.swap) go to 120 c c perform an interchange. c call sswap(imax,a(1,imax),1,a(1,k),1) do 110 jj = imax, k j = k + imax - jj t = a(j,k) a(j,k) = a(imax,j) a(imax,j) = t 110 continue 120 continue c c perform the elimination. c do 130 jj = 1, km1 j = k - jj mulk = -a(j,k)/a(k,k) t = mulk call saxpy(j,t,a(1,k),1,a(1,j),1) a(j,k) = mulk 130 continue c c set the pivot array. c kpvt(k) = k if (swap) kpvt(k) = imax go to 190 140 continue c c 2 x 2 pivot block. c if (.not.swap) go to 160 c c perform an interchange. c call sswap(imax,a(1,imax),1,a(1,k-1),1) do 150 jj = imax, km1 j = km1 + imax - jj t = a(j,k-1) a(j,k-1) = a(imax,j) a(imax,j) = t 150 continue t = a(k-1,k) a(k-1,k) = a(imax,k) a(imax,k) = t 160 continue c c perform the elimination. c km2 = k - 2 if (km2 .eq. 0) go to 180 ak = a(k,k)/a(k-1,k) akm1 = a(k-1,k-1)/a(k-1,k) denom = 1.0e0 - ak*akm1 do 170 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 saxpy(j,t,a(1,k),1,a(1,j),1) t = mulkm1 call saxpy(j,t,a(1,k-1),1,a(1,j),1) a(j,k) = mulk a(j,k-1) = mulkm1 170 continue 180 continue c c set the pivot array. c kpvt(k) = 1 - k if (swap) kpvt(k) = -imax kpvt(k-1) = kpvt(k) 190 continue k = k - kstep go to 10 200 continue return end subroutine ssisl(a,lda,n,kpvt,b) c*********************************************************************72 integer lda,n,kpvt(1) real a(lda,1),b(1) c c ssisl solves the real symmetric system c a * x = b c using the factors computed by ssifa. c c on entry c c a real(lda,n) c the output from ssifa. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c kpvt integer(n) c the pivot vector from ssifa. c c b real(n) c the right hand side vector. c c on return c c b the solution vector x . c c error condition c c a division by zero may occur if ssico has set rcond .eq. 0.0 c or ssifa has set info .ne. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call ssifa(a,lda,n,kpvt,info) c if (info .ne. 0) go to ... c do 10 j = 1, p c call ssisl(a,lda,n,kpvt,c(1,j)) c 10 continue c c linpack. this version dated 08/14/78 . c james bunch, univ. calif. san diego, argonne nat. lab. c c subroutines and functions c c blas saxpy,sdot c fortran iabs c c internal variables. c real ak,akm1,bk,bkm1,sdot,denom,temp integer k,kp c c loop backward applying the transformations and c d inverse to b. c k = n 10 if (k .eq. 0) go to 80 if (kpvt(k) .lt. 0) go to 40 c c 1 x 1 pivot block. c if (k .eq. 1) go to 30 kp = kpvt(k) if (kp .eq. k) go to 20 c c interchange. c temp = b(k) b(k) = b(kp) b(kp) = temp 20 continue c c apply the transformation. c call saxpy(k-1,b(k),a(1,k),1,b(1),1) 30 continue c c apply d inverse. c b(k) = b(k)/a(k,k) k = k - 1 go to 70 40 continue c c 2 x 2 pivot block. c if (k .eq. 2) go to 60 kp = iabs(kpvt(k)) if (kp .eq. k - 1) go to 50 c c interchange. c temp = b(k-1) b(k-1) = b(kp) b(kp) = temp 50 continue c c apply the transformation. c call saxpy(k-2,b(k),a(1,k),1,b(1),1) call saxpy(k-2,b(k-1),a(1,k-1),1,b(1),1) 60 continue c c apply d inverse. c 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 70 continue go to 10 80 continue c c loop forward applying the transformations. c k = 1 90 if (k .gt. n) go to 160 if (kpvt(k) .lt. 0) go to 120 c c 1 x 1 pivot block. c if (k .eq. 1) go to 110 c c apply the transformation. c b(k) = b(k) + sdot(k-1,a(1,k),1,b(1),1) kp = kpvt(k) if (kp .eq. k) go to 100 c c interchange. c temp = b(k) b(k) = b(kp) b(kp) = temp 100 continue 110 continue k = k + 1 go to 150 120 continue c c 2 x 2 pivot block. c if (k .eq. 1) go to 140 c c apply the transformation. c b(k) = b(k) + sdot(k-1,a(1,k),1,b(1),1) b(k+1) = b(k+1) + sdot(k-1,a(1,k+1),1,b(1),1) kp = iabs(kpvt(k)) if (kp .eq. k) go to 130 c c interchange. c temp = b(k) b(k) = b(kp) b(kp) = temp 130 continue 140 continue k = k + 2 150 continue go to 90 160 continue return end subroutine sspco(ap,n,kpvt,rcond,z) c*********************************************************************72 integer n,kpvt(1) real ap(1),z(1) real rcond c c sspco factors a real symmetric matrix stored in packed c form by elimination with symmetric pivoting and estimates c the condition of the matrix. c c if rcond is not needed, sspfa is slightly faster. c to solve a*x = b , follow sspco by sspsl. c to compute inverse(a)*c , follow sspco by sspsl. c to compute inverse(a) , follow sspco by sspdi. c to compute determinant(a) , follow sspco by sspdi. c to compute inertia(a), follow sspco by sspdi. c c on entry c c ap real (n*(n+1)/2) c the packed form of a symmetric matrix a . the c columns of the upper triangle are stored sequentially c in a one-dimensional array of length n*(n+1)/2 . c see comments below for details. c c n integer c the order of the matrix a . c c output c c ap a block diagonal matrix and the multipliers which c were used to obtain it stored in packed form. c the factorization can be written a = u*d*trans(u) c where u is a product of permutation and unit c upper triangular matrices , trans(u) is the c transpose of u , and d is block diagonal c with 1 by 1 and 2 by 2 blocks. c c kpvt integer(n) c an integer vector of pivot indices. c c rcond real c an estimate of the reciprocal condition of a . c for the system a*x = b , relative perturbations c in a and b of size epsilon may cause c relative perturbations in x of size epsilon/rcond . c if rcond is so small that the logical expression c 1.0 + rcond .eq. 1.0 c is true, then a may be singular to working c precision. in particular, rcond is zero if c exact singularity is detected or the estimate c underflows. c c z real(n) c a work vector whose contents are usually unimportant. c if a is close to a singular matrix, then z is c an approximate null vector in the sense that c norm(a*z) = rcond*norm(a)*norm(z) . c c packed storage c c the following program segment will pack the upper c triangle of a symmetric matrix. c c k = 0 c do 20 j = 1, n c do 10 i = 1, j c k = k + 1 c ap(k) = a(i,j) c 10 continue c 20 continue c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c linpack sspfa c blas saxpy,sdot,sscal,sasum c fortran abs,amax1,iabs,sign c c internal variables c real ak,akm1,bk,bkm1,sdot,denom,ek,t real anorm,s,sasum,ynorm integer i,ij,ik,ikm1,ikp1,info,j,jm1,j1 integer k,kk,km1k,km1km1,kp,kps,ks c c c find norm of a using only upper half c j1 = 1 do 30 j = 1, n z(j) = sasum(j,ap(j1),1) ij = j1 j1 = j1 + j jm1 = j - 1 if (jm1 .lt. 1) go to 20 do 10 i = 1, jm1 z(i) = z(i) + abs(ap(ij)) ij = ij + 1 10 continue 20 continue 30 continue anorm = 0.0e0 do 40 j = 1, n anorm = amax1(anorm,z(j)) 40 continue c c factor c call sspfa(ap,n,kpvt,info) c c rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) . c estimate = norm(z)/norm(y) where a*z = y and a*y = e . c the components of e are chosen to cause maximum local c growth in the elements of w where u*d*w = e . c the vectors are frequently rescaled to avoid overflow. c c solve u*d*w = e c ek = 1.0e0 do 50 j = 1, n z(j) = 0.0e0 50 continue k = n ik = (n*(n - 1))/2 60 if (k .eq. 0) go to 120 kk = ik + k ikm1 = ik - (k - 1) ks = 1 if (kpvt(k) .lt. 0) ks = 2 kp = iabs(kpvt(k)) kps = k + 1 - ks if (kp .eq. kps) go to 70 t = z(kps) z(kps) = z(kp) z(kp) = t 70 continue if (z(k) .ne. 0.0e0) ek = sign(ek,z(k)) z(k) = z(k) + ek call saxpy(k-ks,z(k),ap(ik+1),1,z(1),1) if (ks .eq. 1) go to 80 if (z(k-1) .ne. 0.0e0) ek = sign(ek,z(k-1)) z(k-1) = z(k-1) + ek call saxpy(k-ks,z(k-1),ap(ikm1+1),1,z(1),1) 80 continue if (ks .eq. 2) go to 100 if (abs(z(k)) .le. abs(ap(kk))) go to 90 s = abs(ap(kk))/abs(z(k)) call sscal(n,s,z,1) ek = s*ek 90 continue if (ap(kk) .ne. 0.0e0) z(k) = z(k)/ap(kk) if (ap(kk) .eq. 0.0e0) z(k) = 1.0e0 go to 110 100 continue 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 110 continue k = k - ks ik = ik - k if (ks .eq. 2) ik = ik - (k + 1) go to 60 120 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) c c solve trans(u)*y = w c k = 1 ik = 0 130 if (k .gt. n) go to 160 ks = 1 if (kpvt(k) .lt. 0) ks = 2 if (k .eq. 1) go to 150 z(k) = z(k) + sdot(k-1,ap(ik+1),1,z(1),1) ikp1 = ik + k if (ks .eq. 2) * z(k+1) = z(k+1) + sdot(k-1,ap(ikp1+1),1,z(1),1) kp = iabs(kpvt(k)) if (kp .eq. k) go to 140 t = z(k) z(k) = z(kp) z(kp) = t 140 continue 150 continue ik = ik + k if (ks .eq. 2) ik = ik + (k + 1) k = k + ks go to 130 160 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) c ynorm = 1.0e0 c c solve u*d*v = y c k = n ik = n*(n - 1)/2 170 if (k .eq. 0) go to 230 kk = ik + k ikm1 = ik - (k - 1) ks = 1 if (kpvt(k) .lt. 0) ks = 2 if (k .eq. ks) go to 190 kp = iabs(kpvt(k)) kps = k + 1 - ks if (kp .eq. kps) go to 180 t = z(kps) z(kps) = z(kp) z(kp) = t 180 continue call saxpy(k-ks,z(k),ap(ik+1),1,z(1),1) if (ks .eq. 2) call saxpy(k-ks,z(k-1),ap(ikm1+1),1,z(1),1) 190 continue if (ks .eq. 2) go to 210 if (abs(z(k)) .le. abs(ap(kk))) go to 200 s = abs(ap(kk))/abs(z(k)) call sscal(n,s,z,1) ynorm = s*ynorm 200 continue if (ap(kk) .ne. 0.0e0) z(k) = z(k)/ap(kk) if (ap(kk) .eq. 0.0e0) z(k) = 1.0e0 go to 220 210 continue 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 220 continue k = k - ks ik = ik - k if (ks .eq. 2) ik = ik - (k + 1) go to 170 230 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) ynorm = s*ynorm c c solve trans(u)*z = v c k = 1 ik = 0 240 if (k .gt. n) go to 270 ks = 1 if (kpvt(k) .lt. 0) ks = 2 if (k .eq. 1) go to 260 z(k) = z(k) + sdot(k-1,ap(ik+1),1,z(1),1) ikp1 = ik + k if (ks .eq. 2) * z(k+1) = z(k+1) + sdot(k-1,ap(ikp1+1),1,z(1),1) kp = iabs(kpvt(k)) if (kp .eq. k) go to 250 t = z(k) z(k) = z(kp) z(kp) = t 250 continue 260 continue ik = ik + k if (ks .eq. 2) ik = ik + (k + 1) k = k + ks go to 240 270 continue c make znorm = 1.0 s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) ynorm = s*ynorm c if (anorm .ne. 0.0e0) rcond = ynorm/anorm if (anorm .eq. 0.0e0) rcond = 0.0e0 return end subroutine sspdi(ap,n,kpvt,det,inert,work,job) c*********************************************************************72 integer n,job real ap(1),work(1) real det(2) integer kpvt(1),inert(3) c c sspdi computes the determinant, inertia and inverse c of a real symmetric matrix using the factors from sspfa, c where the matrix is stored in packed form. c c on entry c c ap real (n*(n+1)/2) c the output from sspfa. c c n integer c the order of the matrix a. c c kpvt integer(n) c the pivot vector from sspfa. c c work real(n) c work vector. contents ignored. c c job integer c job has the decimal expansion abc where c if c .ne. 0, the inverse is computed, c if b .ne. 0, the determinant is computed, c if a .ne. 0, the inertia is computed. c c for example, job = 111 gives all three. c c on return c c variables not requested by job are not used. c c ap contains the upper triangle of the inverse of c the original matrix, stored in packed form. c the columns of the upper triangle are stored c sequentially in a one-dimensional array. c c det real(2) c determinant of original matrix. c determinant = det(1) * 10.0**det(2) c with 1.0 .le. abs(det(1)) .lt. 10.0 c or det(1) = 0.0. c c inert integer(3) c the inertia of the original matrix. c inert(1) = number of positive eigenvalues. c inert(2) = number of negative eigenvalues. c inert(3) = number of zero eigenvalues. c c error condition c c a division by zero will occur if the inverse is requested c and sspco has set rcond .eq. 0.0 c or sspfa has set info .ne. 0 . c c linpack. this version dated 08/14/78 . c james bunch, univ. calif. san diego, argonne nat. lab. c c subroutines and functions c c blas saxpy,scopy,sdot,sswap c fortran abs,iabs,mod c c internal variables. c real akkp1,sdot,temp real ten,d,t,ak,akp1 integer ij,ik,ikp1,iks,j,jb,jk,jkp1 integer k,kk,kkp1,km1,ks,ksj,kskp1,kstep logical noinv,nodet,noert c noinv = mod(job,10) .eq. 0 nodet = mod(job,100)/10 .eq. 0 noert = mod(job,1000)/100 .eq. 0 c if (nodet .and. noert) go to 140 if (noert) go to 10 inert(1) = 0 inert(2) = 0 inert(3) = 0 10 continue if (nodet) go to 20 det(1) = 1.0e0 det(2) = 0.0e0 ten = 10.0e0 20 continue t = 0.0e0 ik = 0 do 130 k = 1, n kk = ik + k d = ap(kk) c c check if 1 by 1 c if (kpvt(k) .gt. 0) go to 50 c c 2 by 2 block c use det (d s) = (d/t * c - t) * t , t = abs(s) c (s c) c to avoid underflow/overflow troubles. c take two passes through scaling. use t for flag. c if (t .ne. 0.0e0) go to 30 ikp1 = ik + k kkp1 = ikp1 + k t = abs(ap(kkp1)) d = (d/t)*ap(kkp1+1) - t go to 40 30 continue d = t t = 0.0e0 40 continue 50 continue c if (noert) go to 60 if (d .gt. 0.0e0) inert(1) = inert(1) + 1 if (d .lt. 0.0e0) inert(2) = inert(2) + 1 if (d .eq. 0.0e0) inert(3) = inert(3) + 1 60 continue c if (nodet) go to 120 det(1) = d*det(1) if (det(1) .eq. 0.0e0) go to 110 70 if (abs(det(1)) .ge. 1.0e0) go to 80 det(1) = ten*det(1) det(2) = det(2) - 1.0e0 go to 70 80 continue 90 if (abs(det(1)) .lt. ten) go to 100 det(1) = det(1)/ten det(2) = det(2) + 1.0e0 go to 90 100 continue 110 continue 120 continue ik = ik + k 130 continue 140 continue c c compute inverse(a) c if (noinv) go to 270 k = 1 ik = 0 150 if (k .gt. n) go to 260 km1 = k - 1 kk = ik + k ikp1 = ik + k kkp1 = ikp1 + k if (kpvt(k) .lt. 0) go to 180 c c 1 by 1 c ap(kk) = 1.0e0/ap(kk) if (km1 .lt. 1) go to 170 call scopy(km1,ap(ik+1),1,work,1) ij = 0 do 160 j = 1, km1 jk = ik + j ap(jk) = sdot(j,ap(ij+1),1,work,1) call saxpy(j-1,work(j),ap(ij+1),1,ap(ik+1),1) ij = ij + j 160 continue ap(kk) = ap(kk) + sdot(km1,work,1,ap(ik+1),1) 170 continue kstep = 1 go to 220 180 continue c c 2 by 2 c t = abs(ap(kkp1)) ak = ap(kk)/t akp1 = ap(kkp1+1)/t akkp1 = ap(kkp1)/t d = t*(ak*akp1 - 1.0e0) ap(kk) = akp1/d ap(kkp1+1) = ak/d ap(kkp1) = -akkp1/d if (km1 .lt. 1) go to 210 call scopy(km1,ap(ikp1+1),1,work,1) ij = 0 do 190 j = 1, km1 jkp1 = ikp1 + j ap(jkp1) = sdot(j,ap(ij+1),1,work,1) call saxpy(j-1,work(j),ap(ij+1),1,ap(ikp1+1),1) ij = ij + j 190 continue ap(kkp1+1) = ap(kkp1+1) * + sdot(km1,work,1,ap(ikp1+1),1) ap(kkp1) = ap(kkp1) * + sdot(km1,ap(ik+1),1,ap(ikp1+1),1) call scopy(km1,ap(ik+1),1,work,1) ij = 0 do 200 j = 1, km1 jk = ik + j ap(jk) = sdot(j,ap(ij+1),1,work,1) call saxpy(j-1,work(j),ap(ij+1),1,ap(ik+1),1) ij = ij + j 200 continue ap(kk) = ap(kk) + sdot(km1,work,1,ap(ik+1),1) 210 continue kstep = 2 220 continue c c swap c ks = iabs(kpvt(k)) if (ks .eq. k) go to 250 iks = (ks*(ks - 1))/2 call sswap(ks,ap(iks+1),1,ap(ik+1),1) ksj = ik + ks do 230 jb = ks, k j = k + ks - jb jk = ik + j temp = ap(jk) ap(jk) = ap(ksj) ap(ksj) = temp ksj = ksj - (j - 1) 230 continue if (kstep .eq. 1) go to 240 kskp1 = ikp1 + ks temp = ap(kskp1) ap(kskp1) = ap(kkp1) ap(kkp1) = temp 240 continue 250 continue ik = ik + k if (kstep .eq. 2) ik = ik + k + 1 k = k + kstep go to 150 260 continue 270 continue return end subroutine sspfa(ap,n,kpvt,info) c*********************************************************************72 integer n,kpvt(1),info real ap(1) c c sspfa factors a real symmetric matrix stored in c packed form by elimination with symmetric pivoting. c c to solve a*x = b , follow sspfa by sspsl. c to compute inverse(a)*c , follow sspfa by sspsl. c to compute determinant(a) , follow sspfa by sspdi. c to compute inertia(a) , follow sspfa by sspdi. c to compute inverse(a) , follow sspfa by sspdi. c c on entry c c ap real (n*(n+1)/2) c the packed form of a symmetric matrix a . the c columns of the upper triangle are stored sequentially c in a one-dimensional array of length n*(n+1)/2 . c see comments below for details. c c n integer c the order of the matrix a . c c output c c ap a block diagonal matrix and the multipliers which c were used to obtain it stored in packed form. c the factorization can be written a = u*d*trans(u) c where u is a product of permutation and unit c upper triangular matrices , trans(u) is the c transpose of u , and d is block diagonal c with 1 by 1 and 2 by 2 blocks. c c kpvt integer(n) c an integer vector of pivot indices. c c info integer c = 0 normal value. c = k if the k-th pivot block is singular. this is c not an error condition for this subroutine, c but it does indicate that sspsl or sspdi may c divide by zero if called. c c packed storage c c the following program segment will pack the upper c triangle of a symmetric matrix. c c k = 0 c do 20 j = 1, n c do 10 i = 1, j c k = k + 1 c ap(k) = a(i,j) c 10 continue c 20 continue c c linpack. this version dated 08/14/78 . c james bunch, univ. calif. san diego, argonne nat. lab. c c subroutines and functions c c blas saxpy,sswap,isamax c fortran abs,amax1,sqrt c c internal variables c real ak,akm1,bk,bkm1,denom,mulk,mulkm1,t real absakk,alpha,colmax,rowmax integer isamax,ij,ijj,ik,ikm1,im,imax,imaxp1,imim,imj,imk integer j,jj,jk,jkm1,jmax,jmim,k,kk,km1,km1k,km1km1,km2,kstep logical swap c c c initialize c c alpha is used in choosing pivot block size. alpha = (1.0e0 + sqrt(17.0e0))/8.0e0 c info = 0 c c main loop on k, which goes from n to 1. c k = n ik = (n*(n - 1))/2 10 continue c c leave the loop if k=0 or k=1. c c ...exit if (k .eq. 0) go to 200 if (k .gt. 1) go to 20 kpvt(1) = 1 if (ap(1) .eq. 0.0e0) info = 1 c ......exit go to 200 20 continue c c this section of code determines the kind of c elimination to be performed. when it is completed, c kstep will be set to the size of the pivot block, and c swap will be set to .true. if an interchange is c required. c km1 = k - 1 kk = ik + k absakk = abs(ap(kk)) c c determine the largest off-diagonal element in c column k. c imax = isamax(k-1,ap(ik+1),1) imk = ik + imax colmax = abs(ap(imk)) if (absakk .lt. alpha*colmax) go to 30 kstep = 1 swap = .false. go to 90 30 continue c c determine the largest off-diagonal element in c row imax. c rowmax = 0.0e0 imaxp1 = imax + 1 im = imax*(imax - 1)/2 imj = im + 2*imax do 40 j = imaxp1, k rowmax = amax1(rowmax,abs(ap(imj))) imj = imj + j 40 continue if (imax .eq. 1) go to 50 jmax = isamax(imax-1,ap(im+1),1) jmim = jmax + im rowmax = amax1(rowmax,abs(ap(jmim))) 50 continue imim = imax + im if (abs(ap(imim)) .lt. alpha*rowmax) go to 60 kstep = 1 swap = .true. go to 80 60 continue if (absakk .lt. alpha*colmax*(colmax/rowmax)) go to 70 kstep = 1 swap = .false. go to 80 70 continue kstep = 2 swap = imax .ne. km1 80 continue 90 continue if (amax1(absakk,colmax) .ne. 0.0e0) go to 100 c c column k is zero. set info and iterate the loop. c kpvt(k) = k info = k go to 190 100 continue if (kstep .eq. 2) go to 140 c c 1 x 1 pivot block. c if (.not.swap) go to 120 c c perform an interchange. c call sswap(imax,ap(im+1),1,ap(ik+1),1) imj = ik + imax do 110 jj = imax, k j = k + imax - jj jk = ik + j t = ap(jk) ap(jk) = ap(imj) ap(imj) = t imj = imj - (j - 1) 110 continue 120 continue c c perform the elimination. c ij = ik - (k - 1) do 130 jj = 1, km1 j = k - jj jk = ik + j mulk = -ap(jk)/ap(kk) t = mulk call saxpy(j,t,ap(ik+1),1,ap(ij+1),1) ijj = ij + j ap(jk) = mulk ij = ij - (j - 1) 130 continue c c set the pivot array. c kpvt(k) = k if (swap) kpvt(k) = imax go to 190 140 continue c c 2 x 2 pivot block. c km1k = ik + k - 1 ikm1 = ik - (k - 1) if (.not.swap) go to 160 c c perform an interchange. c call sswap(imax,ap(im+1),1,ap(ikm1+1),1) imj = ikm1 + imax do 150 jj = imax, km1 j = km1 + imax - jj jkm1 = ikm1 + j t = ap(jkm1) ap(jkm1) = ap(imj) ap(imj) = t imj = imj - (j - 1) 150 continue t = ap(km1k) ap(km1k) = ap(imk) ap(imk) = t 160 continue c c perform the elimination. c km2 = k - 2 if (km2 .eq. 0) go to 180 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 170 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 saxpy(j,t,ap(ik+1),1,ap(ij+1),1) t = mulkm1 call saxpy(j,t,ap(ikm1+1),1,ap(ij+1),1) ap(jk) = mulk ap(jkm1) = mulkm1 ijj = ij + j ij = ij - (j - 1) 170 continue 180 continue c c set the pivot array. c kpvt(k) = 1 - k if (swap) kpvt(k) = -imax kpvt(k-1) = kpvt(k) 190 continue ik = ik - (k - 1) if (kstep .eq. 2) ik = ik - (k - 2) k = k - kstep go to 10 200 continue return end subroutine sspsl(ap,n,kpvt,b) c*********************************************************************72 integer n,kpvt(1) real ap(1),b(1) c c ssisl solves the real symmetric system c a * x = b c using the factors computed by sspfa. c c on entry c c ap real(n*(n+1)/2) c the output from sspfa. c c n integer c the order of the matrix a . c c kpvt integer(n) c the pivot vector from sspfa. c c b real(n) c the right hand side vector. c c on return c c b the solution vector x . c c error condition c c a division by zero may occur if sspco has set rcond .eq. 0.0 c or sspfa has set info .ne. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call sspfa(ap,n,kpvt,info) c if (info .ne. 0) go to ... c do 10 j = 1, p c call sspsl(ap,n,kpvt,c(1,j)) c 10 continue c c linpack. this version dated 08/14/78 . c james bunch, univ. calif. san diego, argonne nat. lab. c c subroutines and functions c c blas saxpy,sdot c fortran iabs c c internal variables. c real ak,akm1,bk,bkm1,sdot,denom,temp integer ik,ikm1,ikp1,k,kk,km1k,km1km1,kp c c loop backward applying the transformations and c d inverse to b. c k = n ik = (n*(n - 1))/2 10 if (k .eq. 0) go to 80 kk = ik + k if (kpvt(k) .lt. 0) go to 40 c c 1 x 1 pivot block. c if (k .eq. 1) go to 30 kp = kpvt(k) if (kp .eq. k) go to 20 c c interchange. c temp = b(k) b(k) = b(kp) b(kp) = temp 20 continue c c apply the transformation. c call saxpy(k-1,b(k),ap(ik+1),1,b(1),1) 30 continue c c apply d inverse. c b(k) = b(k)/ap(kk) k = k - 1 ik = ik - k go to 70 40 continue c c 2 x 2 pivot block. c ikm1 = ik - (k - 1) if (k .eq. 2) go to 60 kp = iabs(kpvt(k)) if (kp .eq. k - 1) go to 50 c c interchange. c temp = b(k-1) b(k-1) = b(kp) b(kp) = temp 50 continue c c apply the transformation. c call saxpy(k-2,b(k),ap(ik+1),1,b(1),1) call saxpy(k-2,b(k-1),ap(ikm1+1),1,b(1),1) 60 continue c c apply d inverse. c 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 70 continue go to 10 80 continue c c loop forward applying the transformations. c k = 1 ik = 0 90 if (k .gt. n) go to 160 if (kpvt(k) .lt. 0) go to 120 c c 1 x 1 pivot block. c if (k .eq. 1) go to 110 c c apply the transformation. c b(k) = b(k) + sdot(k-1,ap(ik+1),1,b(1),1) kp = kpvt(k) if (kp .eq. k) go to 100 c c interchange. c temp = b(k) b(k) = b(kp) b(kp) = temp 100 continue 110 continue ik = ik + k k = k + 1 go to 150 120 continue c c 2 x 2 pivot block. c if (k .eq. 1) go to 140 c c apply the transformation. c b(k) = b(k) + sdot(k-1,ap(ik+1),1,b(1),1) ikp1 = ik + k b(k+1) = b(k+1) + sdot(k-1,ap(ikp1+1),1,b(1),1) kp = iabs(kpvt(k)) if (kp .eq. k) go to 130 c c interchange. c temp = b(k) b(k) = b(kp) b(kp) = temp 130 continue 140 continue ik = ik + k + k + 1 k = k + 2 150 continue go to 90 160 continue return end subroutine ssvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info) c*********************************************************************72 integer ldx,n,p,ldu,ldv,job,info real x(ldx,1),s(1),e(1),u(ldu,1),v(ldv,1),work(1) c c c ssvdc is a subroutine to reduce a real nxp matrix x by c orthogonal transformations u and v to diagonal form. the c diagonal elements s(i) are the singular values of x. the c columns of u are the corresponding left singular vectors, c and the columns of v the right singular vectors. c c on entry c c x real(ldx,p), where ldx.ge.n. c x contains the matrix whose singular value c decomposition is to be computed. x is c destroyed by ssvdc. c c ldx integer. c ldx is the leading dimension of the array x. c c n integer. c n is the number of rows of the matrix x. c c p integer. c p is the number of columns of the matrix x. c c ldu integer. c ldu is the leading dimension of the array u. c (see below). c c ldv integer. c ldv is the leading dimension of the array v. c (see below). c c work real(n). c work is a scratch array. c c job integer. c job controls the computation of the singular c vectors. it has the decimal expansion ab c with the following meaning c c a.eq.0 do not compute the left singular c vectors. c a.eq.1 return the n left singular vectors c in u. c a.ge.2 return the first min(n,p) singular c vectors in u. c b.eq.0 do not compute the right singular c vectors. c b.eq.1 return the right singular vectors c in v. c c on return c c s real(mm), where mm=max(n+1,p). c the first min(n,p) entries of s contain the c singular values of x arranged in descending c order of magnitude. c c e real(mm), where mm=max(n+1,p). c e ordinarily contains zeros. however see the c discussion of info for exceptions. c c u real(ldu,k), where ldu.ge.n. if joba.eq.1 then c k.eq.n, if joba.ge.2 then c k.eq.min(n,p). c u contains the matrix of left singular vectors. c u is not referenced if joba.eq.0. if n.le.p c or if joba.eq.2, then u may be identified with x c in the subroutine call. c c v real(ldv,p), where ldv.ge.p. c v contains the matrix of right singular vectors. c v is not referenced if job.eq.0. if p.le.n, c then v may be identified with x in the c subroutine call. c c info integer. c the singular values (and their corresponding c singular vectors) s(info+1),s(info+2),...,s(m) c are correct (here m=min(n,p)). thus if c info.eq.0, all the singular values and their c vectors are correct. in any event, the matrix c b = trans(u)*x*v is the bidiagonal matrix c with the elements of s on its diagonal and the c elements of e on its super-diagonal (trans(u) c is the transpose of u). thus the singular c values of x and b are the same. c c linpack. this version dated 03/19/79 . c correction to shift calculation made 2/85. c g.w. stewart, university of maryland, argonne national lab. c c ***** uses the following functions and subprograms. c c external srot c blas saxpy,sdot,sscal,sswap,snrm2,srotg c fortran abs,amax1,max0,min0,mod,sqrt c c internal variables c integer i,iter,j,jobu,k,kase,kk,l,ll,lls,lm1,lp1,ls,lu,m,maxit, * mm,mm1,mp1,nct,nctp1,ncu,nrt,nrtp1 real sdot,t real b,c,cs,el,emm1,f,g,snrm2,scale,shift,sl,sm,sn,smm1,t1,test, * ztest logical wantu,wantv c c c set the maximum number of iterations. c maxit = 30 c c determine what is to be computed. c wantu = .false. wantv = .false. jobu = mod(job,100)/10 ncu = n if (jobu .gt. 1) ncu = min0(n,p) if (jobu .ne. 0) wantu = .true. if (mod(job,10) .ne. 0) wantv = .true. c c reduce x to bidiagonal form, storing the diagonal elements c in s and the super-diagonal elements in e. c info = 0 nct = min0(n-1,p) nrt = max0(0,min0(p-2,n)) lu = max0(nct,nrt) if (lu .lt. 1) go to 170 do 160 l = 1, lu lp1 = l + 1 if (l .gt. nct) go to 20 c c compute the transformation for the l-th column and c place the l-th diagonal in s(l). c s(l) = snrm2(n-l+1,x(l,l),1) if (s(l) .eq. 0.0e0) go to 10 if (x(l,l) .ne. 0.0e0) s(l) = sign(s(l),x(l,l)) call sscal(n-l+1,1.0e0/s(l),x(l,l),1) x(l,l) = 1.0e0 + x(l,l) 10 continue s(l) = -s(l) 20 continue if (p .lt. lp1) go to 50 do 40 j = lp1, p if (l .gt. nct) go to 30 if (s(l) .eq. 0.0e0) go to 30 c c apply the transformation. c t = -sdot(n-l+1,x(l,l),1,x(l,j),1)/x(l,l) call saxpy(n-l+1,t,x(l,l),1,x(l,j),1) 30 continue c c place the l-th row of x into e for the c subsequent calculation of the row transformation. c e(j) = x(l,j) 40 continue 50 continue if (.not.wantu .or. l .gt. nct) go to 70 c c place the transformation in u for subsequent back c multiplication. c do 60 i = l, n u(i,l) = x(i,l) 60 continue 70 continue if (l .gt. nrt) go to 150 c c compute the l-th row transformation and place the c l-th super-diagonal in e(l). c e(l) = snrm2(p-l,e(lp1),1) if (e(l) .eq. 0.0e0) go to 80 if (e(lp1) .ne. 0.0e0) e(l) = sign(e(l),e(lp1)) call sscal(p-l,1.0e0/e(l),e(lp1),1) e(lp1) = 1.0e0 + e(lp1) 80 continue e(l) = -e(l) if (lp1 .gt. n .or. e(l) .eq. 0.0e0) go to 120 c c apply the transformation. c do 90 i = lp1, n work(i) = 0.0e0 90 continue do 100 j = lp1, p call saxpy(n-l,e(j),x(lp1,j),1,work(lp1),1) 100 continue do 110 j = lp1, p call saxpy(n-l,-e(j)/e(lp1),work(lp1),1,x(lp1,j),1) 110 continue 120 continue if (.not.wantv) go to 140 c c place the transformation in v for subsequent c back multiplication. c do 130 i = lp1, p v(i,l) = e(i) 130 continue 140 continue 150 continue 160 continue 170 continue c c set up the final bidiagonal matrix or order m. c m = min0(p,n+1) nctp1 = nct + 1 nrtp1 = nrt + 1 if (nct .lt. p) s(nctp1) = x(nctp1,nctp1) if (n .lt. m) s(m) = 0.0e0 if (nrtp1 .lt. m) e(nrtp1) = x(nrtp1,m) e(m) = 0.0e0 c c if required, generate u. c if (.not.wantu) go to 300 if (ncu .lt. nctp1) go to 200 do 190 j = nctp1, ncu do 180 i = 1, n u(i,j) = 0.0e0 180 continue u(j,j) = 1.0e0 190 continue 200 continue if (nct .lt. 1) go to 290 do 280 ll = 1, nct l = nct - ll + 1 if (s(l) .eq. 0.0e0) go to 250 lp1 = l + 1 if (ncu .lt. lp1) go to 220 do 210 j = lp1, ncu t = -sdot(n-l+1,u(l,l),1,u(l,j),1)/u(l,l) call saxpy(n-l+1,t,u(l,l),1,u(l,j),1) 210 continue 220 continue call sscal(n-l+1,-1.0e0,u(l,l),1) u(l,l) = 1.0e0 + u(l,l) lm1 = l - 1 if (lm1 .lt. 1) go to 240 do 230 i = 1, lm1 u(i,l) = 0.0e0 230 continue 240 continue go to 270 250 continue do 260 i = 1, n u(i,l) = 0.0e0 260 continue u(l,l) = 1.0e0 270 continue 280 continue 290 continue 300 continue c c if it is required, generate v. c if (.not.wantv) go to 350 do 340 ll = 1, p l = p - ll + 1 lp1 = l + 1 if (l .gt. nrt) go to 320 if (e(l) .eq. 0.0e0) go to 320 do 310 j = lp1, p t = -sdot(p-l,v(lp1,l),1,v(lp1,j),1)/v(lp1,l) call saxpy(p-l,t,v(lp1,l),1,v(lp1,j),1) 310 continue 320 continue do 330 i = 1, p v(i,l) = 0.0e0 330 continue v(l,l) = 1.0e0 340 continue 350 continue c c main iteration loop for the singular values. c mm = m iter = 0 360 continue c c quit if all the singular values have been found. c c ...exit if (m .eq. 0) go to 620 c c if too many iterations have been performed, set c flag and return. c if (iter .lt. maxit) go to 370 info = m c ......exit go to 620 370 continue c c this section of the program inspects for c negligible elements in the s and e arrays. on c completion the variables kase and l are set as follows. c c kase = 1 if s(m) and e(l-1) are negligible and l.lt.m c kase = 2 if s(l) is negligible and l.lt.m c kase = 3 if e(l-1) is negligible, l.lt.m, and c s(l), ..., s(m) are not negligible (qr step). c kase = 4 if e(m-1) is negligible (convergence). c do 390 ll = 1, m l = m - ll c ...exit if (l .eq. 0) go to 400 test = abs(s(l)) + abs(s(l+1)) ztest = test + abs(e(l)) if (ztest .ne. test) go to 380 e(l) = 0.0e0 c ......exit go to 400 380 continue 390 continue 400 continue if (l .ne. m - 1) go to 410 kase = 4 go to 480 410 continue lp1 = l + 1 mp1 = m + 1 do 430 lls = lp1, mp1 ls = m - lls + lp1 c ...exit if (ls .eq. l) go to 440 test = 0.0e0 if (ls .ne. m) test = test + abs(e(ls)) if (ls .ne. l + 1) test = test + abs(e(ls-1)) ztest = test + abs(s(ls)) if (ztest .ne. test) go to 420 s(ls) = 0.0e0 c ......exit go to 440 420 continue 430 continue 440 continue if (ls .ne. l) go to 450 kase = 3 go to 470 450 continue if (ls .ne. m) go to 460 kase = 1 go to 470 460 continue kase = 2 l = ls 470 continue 480 continue l = l + 1 c c perform the task indicated by kase. c go to (490,520,540,570), kase c c deflate negligible s(m). c 490 continue mm1 = m - 1 f = e(m-1) e(m-1) = 0.0e0 do 510 kk = l, mm1 k = mm1 - kk + l t1 = s(k) call srotg(t1,f,cs,sn) s(k) = t1 if (k .eq. l) go to 500 f = -sn*e(k-1) e(k-1) = cs*e(k-1) 500 continue if (wantv) call srot(p,v(1,k),1,v(1,m),1,cs,sn) 510 continue go to 610 c c split at negligible s(l). c 520 continue f = e(l-1) e(l-1) = 0.0e0 do 530 k = l, m t1 = s(k) call srotg(t1,f,cs,sn) s(k) = t1 f = -sn*e(k) e(k) = cs*e(k) if (wantu) call srot(n,u(1,k),1,u(1,l-1),1,cs,sn) 530 continue go to 610 c c perform one qr step. c 540 continue c c calculate the shift. c scale = amax1(abs(s(m)),abs(s(m-1)),abs(e(m-1)),abs(s(l)), * abs(e(l))) sm = s(m)/scale smm1 = s(m-1)/scale emm1 = e(m-1)/scale sl = s(l)/scale el = e(l)/scale b = ((smm1 + sm)*(smm1 - sm) + emm1**2)/2.0e0 c = (sm*emm1)**2 shift = 0.0e0 if (b .eq. 0.0e0 .and. c .eq. 0.0e0) go to 550 shift = sqrt(b**2+c) if (b .lt. 0.0e0) shift = -shift shift = c/(b + shift) 550 continue f = (sl + sm)*(sl - sm) + shift g = sl*el c c chase zeros. c mm1 = m - 1 do 560 k = l, mm1 call srotg(f,g,cs,sn) if (k .ne. l) e(k-1) = f f = cs*s(k) + sn*e(k) e(k) = cs*e(k) - sn*s(k) g = sn*s(k+1) s(k+1) = cs*s(k+1) if (wantv) call srot(p,v(1,k),1,v(1,k+1),1,cs,sn) call srotg(f,g,cs,sn) s(k) = f f = cs*e(k) + sn*s(k+1) s(k+1) = -sn*e(k) + cs*s(k+1) g = sn*e(k+1) e(k+1) = cs*e(k+1) if (wantu .and. k .lt. n) * call srot(n,u(1,k),1,u(1,k+1),1,cs,sn) 560 continue e(m-1) = f iter = iter + 1 go to 610 c c convergence. c 570 continue c c make the singular value positive. c if (s(l) .ge. 0.0e0) go to 580 s(l) = -s(l) if (wantv) call sscal(p,-1.0e0,v(1,l),1) 580 continue c c order the singular value. c 590 if (l .eq. mm) go to 600 c ...exit if (s(l) .ge. s(l+1)) go to 600 t = s(l) s(l) = s(l+1) s(l+1) = t if (wantv .and. l .lt. p) * call sswap(p,v(1,l),1,v(1,l+1),1) if (wantu .and. l .lt. n) * call sswap(n,u(1,l),1,u(1,l+1),1) l = l + 1 go to 590 600 continue iter = 0 m = m - 1 610 continue go to 360 620 continue return end subroutine strco(t,ldt,n,rcond,z,job) c*********************************************************************72 integer ldt,n,job real t(ldt,1),z(1) real rcond c c strco estimates the condition of a real triangular matrix. c c on entry c c t real(ldt,n) c t contains the triangular matrix. the zero c elements of the matrix are not referenced, and c the corresponding elements of the array can be c used to store other information. c c ldt integer c ldt is the leading dimension of the array t. c c n integer c n is the order of the system. c c job integer c = 0 t is lower triangular. c = nonzero t is upper triangular. c c on return c c rcond real c an estimate of the reciprocal condition of t . c for the system t*x = b , relative perturbations c in t and b of size epsilon may cause c relative perturbations in x of size epsilon/rcond . c if rcond is so small that the logical expression c 1.0 + rcond .eq. 1.0 c is true, then t may be singular to working c precision. in particular, rcond is zero if c exact singularity is detected or the estimate c underflows. c c z real(n) c a work vector whose contents are usually unimportant. c if t is close to a singular matrix, then z is c an approximate null vector in the sense that c norm(a*z) = rcond*norm(a)*norm(z) . c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas saxpy,sscal,sasum c fortran abs,amax1,sign c c internal variables c real w,wk,wkm,ek real tnorm,ynorm,s,sm,sasum integer i1,j,j1,j2,k,kk,l logical lower c lower = job .eq. 0 c c compute 1-norm of t c tnorm = 0.0e0 do 10 j = 1, n l = j if (lower) l = n + 1 - j i1 = 1 if (lower) i1 = j tnorm = amax1(tnorm,sasum(l,t(i1,j),1)) 10 continue c c rcond = 1/(norm(t)*(estimate of norm(inverse(t)))) . c estimate = norm(z)/norm(y) where t*z = y and trans(t)*y = e . c trans(t) is the transpose of t . c the components of e are chosen to cause maximum local c growth in the elements of y . c the vectors are frequently rescaled to avoid overflow. c c solve trans(t)*y = e c ek = 1.0e0 do 20 j = 1, n z(j) = 0.0e0 20 continue do 100 kk = 1, n k = kk if (lower) k = n + 1 - kk if (z(k) .ne. 0.0e0) ek = sign(ek,-z(k)) if (abs(ek-z(k)) .le. abs(t(k,k))) go to 30 s = abs(t(k,k))/abs(ek-z(k)) call sscal(n,s,z,1) ek = s*ek 30 continue wk = ek - z(k) wkm = -ek - z(k) s = abs(wk) sm = abs(wkm) if (t(k,k) .eq. 0.0e0) go to 40 wk = wk/t(k,k) wkm = wkm/t(k,k) go to 50 40 continue wk = 1.0e0 wkm = 1.0e0 50 continue if (kk .eq. n) go to 90 j1 = k + 1 if (lower) j1 = 1 j2 = n if (lower) j2 = k - 1 do 60 j = j1, j2 sm = sm + abs(z(j)+wkm*t(k,j)) z(j) = z(j) + wk*t(k,j) s = s + abs(z(j)) 60 continue if (s .ge. sm) go to 80 w = wkm - wk wk = wkm do 70 j = j1, j2 z(j) = z(j) + w*t(k,j) 70 continue 80 continue 90 continue z(k) = wk 100 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) c ynorm = 1.0e0 c c solve t*z = y c do 130 kk = 1, n k = n + 1 - kk if (lower) k = kk if (abs(z(k)) .le. abs(t(k,k))) go to 110 s = abs(t(k,k))/abs(z(k)) call sscal(n,s,z,1) ynorm = s*ynorm 110 continue if (t(k,k) .ne. 0.0e0) z(k) = z(k)/t(k,k) if (t(k,k) .eq. 0.0e0) z(k) = 1.0e0 i1 = 1 if (lower) i1 = k + 1 if (kk .ge. n) go to 120 w = -z(k) call saxpy(n-kk,w,t(i1,k),1,z(i1),1) 120 continue 130 continue c make znorm = 1.0 s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) ynorm = s*ynorm c if (tnorm .ne. 0.0e0) rcond = ynorm/tnorm if (tnorm .eq. 0.0e0) rcond = 0.0e0 return end subroutine strdi(t,ldt,n,det,job,info) c*********************************************************************72 integer ldt,n,job,info real t(ldt,1),det(2) c c strdi computes the determinant and inverse of a real c triangular matrix. c c on entry c c t real(ldt,n) c t contains the triangular matrix. the zero c elements of the matrix are not referenced, and c the corresponding elements of the array can be c used to store other information. c c ldt integer c ldt is the leading dimension of the array t. c c n integer c n is the order of the system. c c job integer c = 010 no det, inverse of lower triangular. c = 011 no det, inverse of upper triangular. c = 100 det, no inverse. c = 110 det, inverse of lower triangular. c = 111 det, inverse of upper triangular. c c on return c c t inverse of original matrix if requested. c otherwise unchanged. c c det real(2) c determinant of original matrix if requested. c otherwise not referenced. c determinant = det(1) * 10.0**det(2) c with 1.0 .le. abs(det(1)) .lt. 10.0 c or det(1) .eq. 0.0 . c c info integer c info contains zero if the system is nonsingular c and the inverse is requested. c otherwise info contains the index of c a zero diagonal element of t. c c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas saxpy,sscal c fortran abs,mod c c internal variables c real temp real ten integer i,j,k,kb,km1,kp1 c c begin block permitting ...exits to 180 c c compute determinant c if (job/100 .eq. 0) go to 70 det(1) = 1.0e0 det(2) = 0.0e0 ten = 10.0e0 do 50 i = 1, n det(1) = t(i,i)*det(1) c ...exit if (det(1) .eq. 0.0e0) go to 60 10 if (abs(det(1)) .ge. 1.0e0) go to 20 det(1) = ten*det(1) det(2) = det(2) - 1.0e0 go to 10 20 continue 30 if (abs(det(1)) .lt. ten) go to 40 det(1) = det(1)/ten det(2) = det(2) + 1.0e0 go to 30 40 continue 50 continue 60 continue 70 continue c c compute inverse of upper triangular c if (mod(job/10,10) .eq. 0) go to 170 if (mod(job,10) .eq. 0) go to 120 c begin block permitting ...exits to 110 do 100 k = 1, n info = k c ......exit if (t(k,k) .eq. 0.0e0) go to 110 t(k,k) = 1.0e0/t(k,k) temp = -t(k,k) call sscal(k-1,temp,t(1,k),1) kp1 = k + 1 if (n .lt. kp1) go to 90 do 80 j = kp1, n temp = t(k,j) t(k,j) = 0.0e0 call saxpy(k,temp,t(1,k),1,t(1,j),1) 80 continue 90 continue 100 continue info = 0 110 continue go to 160 120 continue c c compute inverse of lower triangular c do 150 kb = 1, n k = n + 1 - kb info = k c ............exit if (t(k,k) .eq. 0.0e0) go to 180 t(k,k) = 1.0e0/t(k,k) temp = -t(k,k) if (k .ne. n) call sscal(n-k,temp,t(k+1,k),1) km1 = k - 1 if (km1 .lt. 1) go to 140 do 130 j = 1, km1 temp = t(k,j) t(k,j) = 0.0e0 call saxpy(n-k+1,temp,t(k,k),1,t(k,j),1) 130 continue 140 continue 150 continue info = 0 160 continue 170 continue 180 continue return end subroutine strsl(t,ldt,n,b,job,info) c*********************************************************************72 integer ldt,n,job,info real t(ldt,1),b(1) c c c strsl solves systems of the form c c t * x = b c or c trans(t) * x = b c c where t is a triangular matrix of order n. here trans(t) c denotes the transpose of the matrix t. c c on entry c c t real(ldt,n) c t contains the matrix of the system. the zero c elements of the matrix are not referenced, and c the corresponding elements of the array can be c used to store other information. c c ldt integer c ldt is the leading dimension of the array t. c c n integer c n is the order of the system. c c b real(n). c b contains the right hand side of the system. c c job integer c job specifies what kind of system is to be solved. c if job is c c 00 solve t*x=b, t lower triangular, c 01 solve t*x=b, t upper triangular, c 10 solve trans(t)*x=b, t lower triangular, c 11 solve trans(t)*x=b, t upper triangular. c c on return c c b b contains the solution, if info .eq. 0. c otherwise b is unaltered. c c info integer c info contains zero if the system is nonsingular. c otherwise info contains the index of c the first zero diagonal element of t. c c linpack. this version dated 08/14/78 . c g. w. stewart, university of maryland, argonne national lab. c c subroutines and functions c c blas saxpy,sdot c fortran mod c c internal variables c real sdot,temp integer case,j,jj c c begin block permitting ...exits to 150 c c check for zero diagonal elements. c do 10 info = 1, n c ......exit if (t(info,info) .eq. 0.0e0) go to 150 10 continue info = 0 c c determine the task and go to it. c case = 1 if (mod(job,10) .ne. 0) case = 2 if (mod(job,100)/10 .ne. 0) case = case + 2 go to (20,50,80,110), case c c solve t*x=b for t lower triangular c 20 continue b(1) = b(1)/t(1,1) if (n .lt. 2) go to 40 do 30 j = 2, n temp = -b(j-1) call saxpy(n-j+1,temp,t(j,j-1),1,b(j),1) b(j) = b(j)/t(j,j) 30 continue 40 continue go to 140 c c solve t*x=b for t upper triangular. c 50 continue b(n) = b(n)/t(n,n) if (n .lt. 2) go to 70 do 60 jj = 2, n j = n - jj + 1 temp = -b(j+1) call saxpy(j,temp,t(1,j+1),1,b(1),1) b(j) = b(j)/t(j,j) 60 continue 70 continue go to 140 c c solve trans(t)*x=b for t lower triangular. c 80 continue b(n) = b(n)/t(n,n) if (n .lt. 2) go to 100 do 90 jj = 2, n j = n - jj + 1 b(j) = b(j) - sdot(jj-1,t(j+1,j),1,b(j+1),1) b(j) = b(j)/t(j,j) 90 continue 100 continue go to 140 c c solve trans(t)*x=b for t upper triangular. c 110 continue b(1) = b(1)/t(1,1) if (n .lt. 2) go to 130 do 120 j = 2, n b(j) = b(j) - sdot(j-1,t(1,j),1,b(1),1) b(j) = b(j)/t(j,j) 120 continue 130 continue 140 continue 150 continue return end