subroutine cexist(cfile,nrc) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c vah 3:29 am, monday, october 27, 1986 c c query the existence of a file using the inquire statement c nrc=0 if file exists, = -1 if not c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ character*(*) cfile logical lexist integer*4 nrc nrc=-1 inquire(file=cfile,exist=lexist) c print *,' cmscmd,nrc :',cmscmd,nrc if(lexist) nrc=0 return end double precision function chicdf(cs,df,ier) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c constructed from the c c imsl routine name - mdch c c by vah, 16:28:44.31, tue, jul 30, 1991 c c changed into full double precision. calls to mdnor replaced by c calls to pheev c c----------------------------------------------------------------------- c c computer - ibm77: double c c latest revision - 16:30:02.36, tue, jul 30, 1991 c c purpose - chi-squared probability distribution function c c usage - p=chicdf(cs,df,ier) c c arguments cs - input value for which the probability is c computed. cs must be greater than or equal c to zero. c df - input number of degrees of freedom of the c chi-squared distribution. df must be greater c than or equal to .5 and less than or equal c to 200,000. c ier - error parameter. (output) c terminal error c ier = 129 indicates that cs or df was c specified incorrectly. c warning error c ier = 34 indicates that the normal pdf c would have produced an underflow. c c returns: p - output probability that a random variable c which follows the chi-squared distribution c with df degrees of freedom is less than or c equal to cs. c precision/hardware - double/all c c reqd. routines - pheev, dgamma c c----------------------------------------------------------------------- c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c specifications for arguments implicit double precision (a-h,o-z) ccc real cs,df,p double precision cs,df,p c specifications for local variables ccc real pt2 double precision pt2 double precision a,z,dgam,eps,w,w1,b,z1,half,one,thrten,thrd double precision dgamma ccc real rinfm,x,c ccc notice that c is declared but never used, so removed double precision rinfm,x data eps/1.0d-6/,half/5.d-1/,thrten/13.d0/,one/1.d0/ data thrd/.3333333333333333d0/ data pt2/.2222222e0/ ccc data rinfm/-0.7237e+76/ data rinfm/-0.7237d+76/ func(w,a,z)=w*dexp(a*dlog(z)-z) c first executable statement c test for invalid input values if (df .ge. .5 .and. df .le. 2.e5 .and. cs .ge. 0.0) go to 5 ier=129 p=rinfm go to 9000 5 ier=0 c set p=0. if cs is less than or c equal to 10.**(-12) if (cs .gt. 1.e-12) go to 15 10 p=0.0 go to 9005 15 if(df.le.100.) go to 20 c use normal distribution approximation c for large degrees of freedom if(cs.lt.2.0) go to 10 x=((cs/df)**thrd-(one-pt2/df))/sqrt(pt2/df) if (x .gt. 5.0) go to 50 if (x .lt. -18.8055) go to 55 ccc call mdnor (x,p) p=pheev(x) go to 9005 c initialization for calculation using c incomplete gamma function 20 if (cs .gt. 200.) go to 50 a=half*df z=half*cs dgam = dgamma(a) w=dmax1(half*a,thrten) if (z .ge. w) go to 35 if (df .gt. 25. .and. cs .lt. 2.) go to 10 c calculate using equation no. 6.5.29 w=one/(dgam*a) w1=w do 25 i=1,50 b=i w1=w1*z/(a+b) if (w1 .le. eps*w) go to 30 w=w+w1 25 continue 30 p=func(w,a,z) go to 9005 c calculate using equation no. 6.5.32 35 z1=one/z b=a-one w1=b*z1 w=one+w1 do 40 i=2,50 b=b-one w1=w1*b*z1 if (w1 .le. eps*w) go to 45 w=w+w1 40 continue 45 w=z1*func(w,a,z) p=one-w/dgam go to 9005 50 p=1.0 go to 9005 c warning error - underflow would have c occurred 55 p=0.0 ier=34 9000 continue ccc call uertst (ier,'mdch ') ccc 9005 return 9005 continue chicdf=p return end subroutine daxpy(n,da,dx,incx,dy,incy) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c constant times a vector plus a vector. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ double precision dx(1),dy(1),da integer i,incx,incy,ixiy,m,mp1,n c if(n.le.0)return if (da .eq. 0.0d0) return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dy(i) + da*dx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 dy(i) = dy(i) + da*dx(i) dy(i + 1) = dy(i + 1) + da*dx(i + 1) dy(i + 2) = dy(i + 2) + da*dx(i + 2) dy(i + 3) = dy(i + 3) + da*dx(i + 3) 50 continue return end subroutine dchdc(a,lda,p,work,jpvt,job,info) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c integer lda,p,jpvt(1),job,info double precision a(lda,1),work(1) c c dchdc 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 double precision(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 double precision. 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 08/14/78 . c j.j. dongarra and g.w. stewart, argonne national laboratory and c university of maryland. c c c blas daxpy,dswap c fortran dsqrt c c internal variables c integer pu,pl,plp1,i,j,jp,jt,k,kb,km1,kp1,l,maxl double precision temp double precision 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 dswap(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 dswap(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) .gt. 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.0d0) 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 dswap(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) = dsqrt(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 daxpy(j-k,temp,work(kp1),1,a(kp1,j),1) 250 continue 260 continue 270 continue 280 continue return end subroutine dcopy(n,dx,incx,dy,incy) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c copies a vector, x, to a vector, y. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ double precision dx(1),dy(1) integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,7) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dx(i) 30 continue if( n .lt. 7 ) return 40 mp1 = m + 1 do 50 i = mp1,n,7 dy(i) = dx(i) dy(i + 1) = dx(i + 1) dy(i + 2) = dx(i + 2) dy(i + 3) = dx(i + 3) dy(i + 4) = dx(i + 4) dy(i + 5) = dx(i + 5) dy(i + 6) = dx(i + 6) 50 continue return end double precision function ddot(n,dx,incx,dy,incy) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c forms the dot product of two vectors. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ double precision dx(1),dy(1),dtemp integer i,incx,incy,ix,iy,m,mp1,n c ddot = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dtemp + dx(ix)*dy(iy) ix = ix + incx iy = iy + incy 10 continue ddot = dtemp return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dx(i)*dy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,5 dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) 50 continue 60 ddot = dtemp return end function dgamma (x) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c specifications for arguments double precision x,dgamma c specifications for local variables dimension p(9),q(8),p4(7) double precision a,b,den,p,p4,q,pi,t,top,big1, 1 xinf,xmin,y,sign,r integer i,ier,j logical mflag c coefficients for minimax c approximation to gamma(x), c 2.0 .le. x .le. 3.0 data p(1)/-5.966047488753637d01/, * p(2)/5.864023793062003d01/, * p(3)/-1.364106217165365d03/, * p(4)/-8.117569271425580d02/, * p(5)/-1.569414683149179d04/, * p(6)/-1.525979925758372d04/, * p(7)/-7.264059615964330d04/, * p(8)/-8.972275718101010d-01/, * p(9)/3.349618189847578d00/ data q(1)/4.103991474182904d02/, * q(2)/-2.262590291514875d03/, * q(3)/2.494325576714903d03/, * q(4)/2.362106244383048d04/, * q(5)/-5.741873227396418d04/, * q(6)/-7.257239715408240d04/, * q(7)/-9.491399521686949d00/, * q(8)/-3.255006939455704d01/ c coefficients for minimax c approximation to ln(gamma(x)), c 12.0 .le. x data p4(1)/8.40596949829d-04/, 1 p4(2)/-5.9523334141881d-04/, 2 p4(3)/7.9365078409227d-04/, * p4(4)/-2.777777777769526d-03/, * p4(5)/8.333333333333333d-02/, * p4(6)/9.189385332046727d-01/, 6 p4(7)/-1.7816839846d-03/ data pi/3.141592653589793d0/ data xinf/1.7d+38/ c dgamma(xmin) .approx. xinf c dgamma(big1) .approx. xinf data xmin/5.8775d-39/ data big1/34.844d0/ c first executable statement ier = 0 mflag = .false. t = x if (dabs(t).gt.xmin) go to 5 ier = 130 dgamma = xinf if (t.le.0.0d0) dgamma = -xinf go to 9000 5 if (dabs(t).lt.big1) go to 10 ier = 129 dgamma = xinf go to 9000 10 if (t.gt.0.0d0) go to 25 c argument is negative mflag = .true. t = -t r = dint(t) sign = 1.0d0 if (dmod(t,2.0d0).eq.0.0d0) sign = -1.0d0 r = t-r if (r.ne.0.0d0) go to 20 ier = 130 dgamma = xinf if (sign.eq.-1.0d0) dgamma = -xinf go to 9000 c argument is not a negative integer 20 r = pi/dsin(r*pi)*sign t = t+1.0d0 c evaluate approximation for dgamma(t) c t .gt. xmin 25 if (t.gt.12.0d0) go to 60 i = t a = 1.0d0 if (i.gt.2) go to 40 i = i+1 go to (30,35,50),i c 0.0 .lt. t .lt. 1.0 30 a = a/(t*(t+1.0d0)) t = t+2.0d0 go to 50 c 1.0 .le. t .lt. 2.0 35 a = a/t t = t+1.0d0 go to 50 c 3.0 .le. t .le. 12.0 40 do 45 j=3,i t = t-1.0d0 a = a*t 45 continue c 2.0 .le. t .le. 3.0 50 top = p(8)*t+p(9) den = t+q(8) do 55 j=1,7 top = top*t+p(j) den = den*t+q(j) 55 continue y = (top/den)*a if (mflag) y = r/y dgamma = y go to 9005 c t .gt. 12.0 60 top = dlog(t) top = t*(top-1.0d0)-.5d0*top t = 1.0d0/t b = t*t a = p4(7) do 65 j = 1,5 65 a = a*b+p4(j) y = a*t+p4(6)+top y = dexp(y) if (mflag) y = r/y dgamma = y go to 9005 9000 continue c call uertst(-ier,6hmgamad) c call uertst(ier,6hdgamma) 9005 return end subroutine dsidi(a,lda,n,kpvt,det,inert,work,job) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c integer lda,n,job double precision a(lda,1),work(1) double precision det(2) integer kpvt(1),inert(3) c c dsidi computes the determinant, inertia and inverse c of a double precision symmetric matrix using the factors from c dsifa. c c on entry c c a double precision(lda,n) c the output from dsifa. 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 dsifa. c c work double precision(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 double precision(2) c determinant of original matrix. c determinant = det(1) * 10.0**det(2) c with 1.0 .le. dabs(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 dsico has set rcond .eq. 0.0 c or dsifa 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 daxpy,dcopy,ddot,dswap c fortran dabs,iabs,mod c c internal variables. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ double precision akkp1,ddot,temp double precision 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.0d0 det(2) = 0.0d0 ten = 10.0d0 20 continue t = 0.0d0 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 = dabs(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.0d0) go to 30 t = dabs(a(k,k+1)) d = (d/t)*a(k+1,k+1) - t go to 40 30 continue d = t t = 0.0d0 40 continue 50 continue c if (noert) go to 60 if (d .gt. 0.0d0) inert(1) = inert(1) + 1 if (d .lt. 0.0d0) inert(2) = inert(2) + 1 if (d .eq. 0.0d0) inert(3) = inert(3) + 1 60 continue c if (nodet) go to 120 det(1) = d*det(1) if (det(1) .eq. 0.0d0) go to 110 70 if (dabs(det(1)) .ge. 1.0d0) go to 80 det(1) = ten*det(1) det(2) = det(2) - 1.0d0 go to 70 80 continue 90 if (dabs(det(1)) .lt. ten) go to 100 det(1) = det(1)/ten det(2) = det(2) + 1.0d0 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.0d0/a(k,k) if (km1 .lt. 1) go to 170 call dcopy(km1,a(1,k),1,work,1) do 160 j = 1, km1 a(j,k) = ddot(j,a(1,j),1,work,1) call daxpy(j-1,work(j),a(1,j),1,a(1,k),1) 160 continue a(k,k) = a(k,k) + ddot(km1,work,1,a(1,k),1) 170 continue kstep = 1 go to 220 180 continue c c 2 by 2 c t = dabs(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.0d0) 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 dcopy(km1,a(1,k+1),1,work,1) do 190 j = 1, km1 a(j,k+1) = ddot(j,a(1,j),1,work,1) call daxpy(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) + ddot(km1,work,1,a(1,k+1),1) a(k,k+1) = a(k,k+1) + ddot(km1,a(1,k),1,a(1,k+1),1) call dcopy(km1,a(1,k),1,work,1) do 200 j = 1, km1 a(j,k) = ddot(j,a(1,j),1,work,1) call daxpy(j-1,work(j),a(1,j),1,a(1,k),1) 200 continue a(k,k) = a(k,k) + ddot(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 dswap(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 dsifa(a,lda,n,kpvt,info) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c integer lda,n,kpvt(1),info double precision a(lda,1) c c dsifa factors a double precision symmetric matrix by elimination c with symmetric pivoting. c c to solve a*x = b , follow dsifa by dsisl. c to compute inverse(a)*c , follow dsifa by dsisl. c to compute determinant(a) , follow dsifa by dsidi. c to compute inertia(a) , follow dsifa by dsidi. c to compute inverse(a) , follow dsifa by dsidi. c c on entry c c a double precision(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 dsisl or dsidi 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 daxpy,dswap,idamax c fortran dabs,dmax1,dsqrt c c internal variables c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ double precision ak,akm1,bk,bkm1,denom,mulk,mulkm1,t double precision absakk,alpha,colmax,rowmax integer imax,imaxp1,j,jj,jmax,k,km1,km2,kstep,idamax logical swap c c c initialize c c alpha is used in choosing pivot block size. alpha = (1.0d0 + dsqrt(17.0d0))/8.0d0 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.0d0) 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 = dabs(a(k,k)) c c determine the largest off-diagonal element in c column k. c imax = idamax(k-1,a(1,k),1) colmax = dabs(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.0d0 imaxp1 = imax + 1 do 40 j = imaxp1, k rowmax = dmax1(rowmax,dabs(a(imax,j))) 40 continue if (imax .eq. 1) go to 50 jmax = idamax(imax-1,a(1,imax),1) rowmax = dmax1(rowmax,dabs(a(jmax,imax))) 50 continue if (dabs(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 (dmax1(absakk,colmax) .ne. 0.0d0) 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 dswap(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 daxpy(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 dswap(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.0d0 - 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 daxpy(j,t,a(1,k),1,a(1,j),1) t = mulkm1 call daxpy(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 dswap (n,dx,incx,dy,incy) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c interchanges two vectors. c uses unrolled loops for increments equal one. c jack dongarra, linpack, 3/11/78. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ double precision dx(1),dy(1),dtemp integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments not equal c to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dx(ix) dx(ix) = dy(iy) dy(iy) = dtemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,3) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp 30 continue if( n .lt. 3 ) return 40 mp1 = m + 1 do 50 i = mp1,n,3 dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp dtemp = dx(i + 1) dx(i + 1) = dy(i + 1) dy(i + 1) = dtemp dtemp = dx(i + 2) dx(i + 2) = dy(i + 2) dy(i + 2) = dtemp 50 continue return end subroutine f3swap(x,y) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c required for pheev3 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision(a-h,o-z) z=x x=y y=z return end subroutine fillm(vec,npi,npj,val) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ double precision vec(npi,npj),val do 100 j=1,npj do 101 i=1,npi vec(i,j)=val 101 continue 100 continue return end subroutine fillv(a,np,val) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ double precision a(np),val do 1 i=1,np a(i)=val 1 continue return end double precision function fmax(vec,n) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision(a-h,o-z) dimension vec(n) tempmax=-1.d20 do 1 iii=1,n temp=vec(iii) if(temp .gt. tempmax) tempmax=temp 1 continue fmax=tempmax return end double precision function fmin(vec,n) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision(a-h,o-z) dimension vec(n) tempmin=1.d20 do 1 iii=1,n temp=vec(iii) if(temp .lt. tempmin) tempmin=temp 1 continue fmin=tempmin return end double precision function hsec() c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision (a-h,o-z) external sec hsec = 100.d0 * sec() return end double precision function hsec1990() c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c vah 17:25:23.43, sat, nov 24, 1990 c c returns number of hundredth's of a second since jan 1, 1990 c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision(a-h,o-z) dimension idaysinm(12),idate(3),itime(5) external hsec,timdat idaysinm(1)=31 idaysinm(2)=28 idaysinm(3)=31 idaysinm(4)=30 idaysinm(5)=31 idaysinm(6)=30 idaysinm(7)=31 idaysinm(8)=31 idaysinm(9)=30 idaysinm(10)=31 idaysinm(11)=30 idaysinm(12)=31 call timdat(itime,idate) idom=idate(1) month=idate(2) year=idate(3) if(year/100.d0 .lt. 1.d0) year=year+1900.d0 if(dmod(year,4.d0) .eq. 0.d0) idaysinm(2)=29 doy=idom do 1 i=1,month-1 doy=doy+idaysinm(i) 1 continue hsec1990=hsec()+100.*3600.*24.*((year-1990.)*365.+doy) return end integer function icoccfst(needle,haystack) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c ----------------------------------------------------------- C c c vah 3/10/1988,10:51:56:23 c 20:40:48.97, Sat, Sep 1, 1990 c c Returns integer location where needle occurs first in haystack. c Returns 0 if no match. c C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ character*(*) needle,haystack lenn=len(needle) lenh=len(haystack) nnn=0 maxlen=lenh-lenn+1 do 1 i=1,maxlen cdbg print *,i,needle(1:lenn),haystack(i:(i+lenn-1)), cdbg * (needle(1:lenn).eq.haystack(i:(i+lenn-1))) if(needle(1:lenn).eq.haystack(i:(i+lenn-1))) then nnn=i goto 2 endif 1 continue 2 continue icoccfst=nnn cdbg print *,'icoccfst returns',nnn,needle,haystack cdbg pause return end integer function icocclst(needle,haystack) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c C ----------------------------------------------------------- C c c vah 14:30:53.47, Sat, Sep 1, 1990 c c Returns integer location where needle occurs last in haystack. c Returns 0 if no match. c C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ character*(*) needle,haystack lenn=len(needle) lenh=len(haystack) nnn=0 maxlen=lenh-lenn+1 do 1 i=1,maxlen cdbg print *,i,needle(1:lenn),haystack(lenh-lenn+2-i:lenh-i+1), cdbg * (needle(1:lenn).eq.haystack(lenh-lenn+2-i:lenh-i+1)) if(needle(1:lenn).eq.haystack(lenh-lenn+2-i:lenh-i+1)) then nnn=lenh-lenn+2-i goto 2 endif 1 continue 2 continue icocclst=nnn cdbg print *,'icocclst returns',nnn,needle,haystack cdbg pause return end integer function idamax(n,dx,incx) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ double precision dx(1),dmax integer i,incx,ix,n c idamax = 0 if( n .lt. 1 ) return idamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 dmax = dabs(dx(1)) ix = ix + incx do 10 i = 2,n if(dabs(dx(ix)).le.dmax) go to 5 idamax = i dmax = dabs(dx(ix)) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 dmax = dabs(dx(1)) do 30 i = 2,n if(dabs(dx(i)).le.dmax) go to 30 idamax = i dmax = dabs(dx(i)) 30 continue return c======================================================================= c dlinpack routines for inversion end c======================================================================= end subroutine invv(a,b,ia,ib,c,n1,npr,nzr,nnr,rcond,*) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c vah : 6:14 pm, friday, july 11, 1986 c*********************************************************************** c c nb: assumes a is symmetric c ********* c c on entry, c ***** c a contains the matrix to be inverted, c ia is the lead dim of a, ib the lead dim of b, c n1 the effective dim. c c on exit, c **** c a is the inverse, b contains the eigenvectors, c c is the vector of eigenvalues in descending order, c npr=# of +ve roots, nnr=# of -ve roots, nzr=# of 0 roots c c error-returns c ************ c case 1: eigen value calculations did not converge(matev2 int=2) c case 2: a root is smaller or equal c to machine absolute precision (eps2), without attempting to c invert the matrix. in both cases, a will be wrong in general c c written by vassilis hajivassiliou 30mar85, c using matev2(a,b,n1,ia,ib,int) c from gqopt. c c eps1= machine relative accuracy (max # of significant digits) c eps2= machine absolute accuracy (smallest epsilon > 0) c*********************************************************************** c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision(a-h,o-z) dimension a(ia,n1),b(ib,n1),c(n1) ccc real*4 x logical novec data eps1/1.d-12/ data eps2/1.d-50/ ccc common/bprec/eps1,eps2 int=1 c since we'll need the eigen vectors novec = int .eq. 0 1 n = n1 xn = n if(n-2) 2,13,3 2 if(.not. novec) b(1,1) = 1.d0 return 3 nc = n - 2 do 12 k = 1,nc kp = k + 1 hq=0.d0 do 9 ml1=kp,n 9 hq=hq+a(ml1,k)*a(ml1,k) s=dsign(dsqrt(hq),a(k+1,k)) if(dabs(s) .lt. 1.d-30) go to 12 a(k+1,k) = a(k+1,k) + s alpha = 1.d0/dsqrt(s * a(k+1,k)) do 4 kr = kp,n 4 a(kr,k) = alpha*a(kr, k) do 5 kr = kp, n hq=0.d0 do 10 ml2=kp,n 10 hq=hq+a(kr,ml2)*a(ml2,k) 5 a(k,kr)=hq hq=0.d0 do 11 ml3=kp,n 11 hq=hq+a(k,ml3)*a(ml3,k) beta=.5d0*hq do 6 kr = kp,n 6 a(k,kr) = a(k,kr) - beta*a(kr,k) do 8 kr = kp,n do 7 j = kr,n a(kr,j) = a(kr, j)-a(k,j)*a(kr,k) - a(k,kr)*a(j,k) 7 a(j,kr) = a(kr, j) 8 continue 12 a(k,k+1) = -s 13 continue hq=0.d0 do 14 ml4=1,n 14 hq=hq+a(ml4,ml4)*a(ml4,ml4) nn=n-1 do 15 ml9=1,nn 15 hq=hq+2.d0*a(ml9,ml9+1)*a(ml9,ml9+1) s=dsqrt(hq) eps = dmax1(s*eps1/dsqrt(xn),eps2) if (novec) go to 280 do 210 k = 1,n do 208 j = 1,n 208 b(j,k) = 0.d0 210 b(k,k) = 1.d0 if (n .eq. 2) go to 280 do 230 k = 1,nc kp = k + 1 do 225 kr = 2,n hq=0.d0 do 228 ml5=kp,n 228 hq=hq+b(kr,ml5)*a(ml5,k) s=hq do 225 j = kp,n 225 b(kr,j) = b(kr,j) - s*a(j,k) 230 continue 280 continue nl = n 304 nf = 1 305 if(nf .ne. nl) go to 308 306 if(nf .eq. 1) go to 999 nl = nf - 1 go to 304 308 nd = nl - nf nda = nl - 1 itcnt = 0 itlim = 30*(nl+1-nf) 309 do 310 kk = 1,nd k = nl - kk 3100 if(dabs(a(k,k+1)).le. eps) go to 315 310 continue go to 322 315 nf = k + 1 go to 305 322 q = a(nf, nf+1) p = a(nf,nf) - 0.5d0*(a(nl,nl) + a(nl-1,nl-1)) - dsign( x dsqrt(.25d0*(a(nl,nl)-a(nl-1,nl-1))**2 + a(nl-1,nl)**2) , y a(nl,nl) + a(nl-1,nl-1) ) itcnt = itcnt + 1 do 350 k = nf,nda if(k .eq. nf) go to 325 p = a(k-1,k) q = r 325 s = p**2 + q**2 pp = p**2/s qq = q**2/s pq = p*q/s t = 2.d0*pq*a(k, k+1) a(k,k+1) = (pp-qq)*a(k,k+1) + pq*(a(k+1,k+1) - a(k,k) ) u = a(k,k) a(k,k) = pp*u + qq*a(k+1,k+1) + t a(k+1,k+1) = qq*u + pp*a(k+1,k+1) - t s = dsqrt(s) p = p/s q = q/s if(k .ne. nf) a(k-1,k) = p*a(k-1,k) + q*r if(k .eq. nda) goto 340 r = q*a(k+1,k+2) a(k+1,k+2) = p*a(k+1,k+2) 340 if(novec) go to 350 do 348 j = 1,n t = p*b(j,k) + q*b(j,k+1) b(j,k+1) =-q*b(j,k) + p*b(j,k+1) 348 b(j,k) = t 350 continue if(itcnt .le. itlim) go to 309 int = 2 go to 306 999 continue nd = n-1 do 800 k = 1,nd s = a(k,k) kr = k kp = kr+1 do 710 j = kp,n if(a(j,j).le.s) go to 710 kr = j s = a(kr,kr) 710 continue if(kr .eq.k) go to 800 a(kr,kr) = a(k,k) a(k,k) = s if (novec) go to 800 do 725 j=1,n t = b(j,k) b(j,k) = b(j,kr) 725 b(j,kr) = t 800 continue c check convergence of roots: if(int.eq.2) then write(*,*) ' eigen values calculations did not converge' write(*,*) ' inverse maybe unreliable' return 1 endif c copy eigenvalues: do 161 k=1,n1 161 c(k)=a(k,k) npr=0 nzr=0 nnr=0 do 171 i=1,n1 if(c(i).gt.0.d0) then npr=npr+1 elseif(c(i).lt.0.d0) then nnr=nnr+1 else nzr=nzr+1 endif 171 continue rcond=-999.d0 do 160 i=1,n1 do 160 j=1,i a(i,j)=0.d0 do 155 k=1,n1 if(dabs(c(k)).le.eps2) goto 979 155 a(i,j)=a(i,j)+b(i,k)*b(j,k)/c(k) 160 a(j,i)=a(i,j) c now we know that roots are fine: rcond=c(1)/c(n1) return 979 write(*,998) k,c(k) 998 format(' at least one root smaller than machine accuracy:'/ * ' root(',i2,')=',d15.7/' no inversion attempted') return 1 end function iseednew(iseed0,icycle) integer*4 iseednew c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c vah 14:38:37.59, Mon, Mar 16, 1992 c c Starts from iseed0, and calls randu icycle times. c For the same iseed0 value, it will thus return a different seed c for every different value of icycle. c c iseed0 and icycle are NOT changed. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ double precision randu,temp integer*4 itemp,iseed0 itemp=iseed0 do 1 i=1,icycle temp=randu(itemp) 1 continue iseednew=itemp return end subroutine lc(ch) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c changes the character variable ch to lower case c vah 3:26 am, monday, october 27, 1986 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ integer*4 length character*(*) ch character*1 c1 length=len(ch) do 1 i=1,length c1=ch(i:i) ic=ichar(c1) c print *,'ic:',ic if(ic.gt.64.and.ic.lt.91) ch(i:i)=char(ic+32) 1 continue return end integer function ltrim(string) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c vah 6:47 pm, saturday, february 27, 1988 character*(*) string c c return the blank-stripped length c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ do 10 l=len(string),1,-1 ltrim=l if (string(l:l) .ne. ' ') return 10 continue return end subroutine machine(machnam) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c vah 8:41:23.91, Tue, Mar 24, 1992 c c returns the value of the environment variable 'machine' c requires machine-dependent subprogram GETEVAR c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ character*(*) machnam call getevar('machine',machnam) return end subroutine matptv(a,lda,b,ldb,c,ldc,n,k,l) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c c=a-transpose*b c c a(lda,k) b(ldb,l) c(ldc,l) c nxk nxl kxl c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ double precision a(lda,k),b(ldb,l),c(ldc,l),hq do 1 i=1,k do 2 j=1,l hq=0.d00 do 3 ijk=1,n hq=hq+a(ijk,i)*b(ijk,j) 3 continue c(i,j)=hq 2 continue 1 continue return end subroutine matpv(a,lda,b,ldb,c,ldc,n,k,l) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c c=a*b c c a(lda,k) b(ldb,l) c(ldc,l) c nxk kxl nxl c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ double precision a(lda,k),b(ldb,l),c(ldc,l),hq do 1 i=1,n do 2 j=1,l hq=0.d00 do 3 ijk=1,k hq=hq+a(i,ijk)*b(ijk,j) 3 continue c(i,j)=hq 2 continue 1 continue return end subroutine mprint(a,lda,n,m,aname,nfile,ndig) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c vah : note that aname is character*(*) c c dbpm trace print of matrix c prints the matrix a in table of n rows and m columns c lda is the no. of rows in dimension of a c aname is the name of array c nfile is the file unit number where the output should go c ndig=6,7 or 8 for 6 per line in d16.8 c ndig=3,4 or 5 for 8 per line in d12.4 c ndig=0,1 or 2 for 10 per line in d9.1 c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision (a-h,o-z) character*(*) aname dimension a(lda,m),outst(10) write(nfile,1) aname,n,m 1 format(/5x,' array ',a8,' (',i4,' by ',i4,' )'/) j1=1 c max no. of columns jmax=max0(10-(ndig/3)*2,6) jj=min0(jmax,m) 100 do 20 i=1,n do 21 j=j1,jj 21 outst(j-j1+1)=a(i,j) if(jmax.ge.10) * write(nfile,10)(outst(j-j1+1),j=j1,jj) if(jmax.ge.8.and.jmax.lt.10) * write(nfile,8)(outst(j-j1+1),j=j1,jj) if(jmax.lt.8) * write(nfile,6)(outst(j-j1+1),j=j1,jj) 6 format(6 (1pd16.8)) 8 format(8 (1pd12.4)) 10 format(10(1pd 9.1)) c 20 continue if(jj.ge.m)write(nfile,2) if(jj.ge.m)return j1=jj+1 jj=min0(jj+jmax,m) if(n.gt.1)write(nfile,2) 2 format(/) goto 100 end double precision function nrmpdf(x,xm,sd) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c ----------------------------------------------------------- c c vah 2:48 pm, friday, april 18, 1986 c c 11:38 pm, wednesday, september 24, 1986 c c msl 4:20 pm, thursday, october 6, 1988 c c vah+msl, 5:40 pm, wed, dec 7, 1988 (pareto pdf) c c ----------------------------------------------------------- c c ----------------------------------------------------------- c c c c 1. double precision function nrmpdf(x,xm,sd) c c 2. double precision function unfpdf(x,a,b) c c 3. double precision function betpdf(arg,a,d) c c 4. double precision function binpdf(x,t,p) c c 5. double precision function caupdf(x,xloc,scale) c c 6. double precision function combnr(xn,xr) c c 7. double precision function exppdf(x,theta) c c 8. double precision function gampdf(x,theta,t) c c 9. double precision function parpdf(x,theta) c c 10. double precision function permnr(xn,xr) c c 11. double precision function dgamcp(x) c c 12. double precision function xfact(x) c c 13. double precision function logpdf(arg,xm,std) c c 14. double precision function dexpdf(arg,xm,std) c c c c ----------------------------------------------------------- c implicit double precision (a-h,o-z) implicit integer*4 (i-n) external sfee c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ arg=(x-xm)/sd nrmpdf=sfee(arg)/sd return end double precision function ph3del(x,y) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c required for pheev3 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision(a-h,o-z) if(x*y)1,3,2 1 ph3del=1. return 2 ph3del=0. return 3 if(x+y)1,2,2 end double precision function ph3s(x,a,b,phix,txa) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c required for pheev3 implicit double precision(a-h,o-z) c computes steck's (1958 ann math stat) s-fn with max error e-9, c though computer rounding errors may reduce this accuracy. calls c function subprograms ph3stk(x,a,b) in which dabs(b) < 1.1, c pheev(x), and ph3t(x,a,phix) note that phix = pheev(x) c and txa = ph3t(x,a,phix). ba = dabs(b) if(ba.gt.1.1)go to 1 ph3s=ph3stk(x,a,b) return 1 ta=dabs(txa) aa=dabs(a) biv=1./ba xf=dabs(phix-.5) xa=dabs(x) xab=xa*aa*ba phixab=pheev(xab) xfa=phixab-.5 if(aa.le.1.)go to 2 ph3s= * xfa*ta-xf*ph3t(xab,biv,phixab)+.25*(xf-xfa)+ph3stk(xab,biv,1./aa) go to 3 2 ab=aa*ba ph3s=.25*(xf+.5)+xfa*ta-ph3stk(xab,1./ab,aa)-ph3stk(xa,ab,biv) 3 if(x.ge.0.)go to 4 ph3s=.1591549*datan(ba/dsqrt(1.+a*a*(1.+b*b)))-s 4 if(b.ge.0.)return ph3s=-ph3s return end double precision function ph3stk(x,a,b) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c required for pheev3 implicit double precision(a-h,o-z) c accurate to .0001 z=1.+a*a z1=dsqrt(z+(a*.5*b)**2) z2=b/dsqrt(z+(a*b)**2) ph3stk=.1591549*pheev(x*z1)*datan(z2) return end double precision function ph3t(z,a,phix) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c required for pheev3 implicit double precision(a-h,o-z) c owens t double precision function accurate to .00005 c phix is a dummy argument x1=z a1=a ff1=0. f1=1. ph3t=0. if(a .ge. 0.)go to 10 a1=-a1 f1=-1. 10 if(a1 .le. 1.)go to 15 x1=dabs(x1) gg1=pheev(x1) x1=a1*x1 a1=1./a1 gg2=pheev(x1) ff1=.5*(gg1+gg2)-gg1*gg2 ff1=ff1*f1 f1=-f1 15 theta=datan(a1) x2a=x1*x1 if(x2a .gt. 150.)go to 20 ph3t=f1*theta*.15915494*dexp(-.5*x2a*a1/theta) ph3t=ph3t*(1.+.00868*(x1*a1)**4) ph3t=ph3t+ff1 20 return end double precision function pheev(y) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c this essentially uses the imsl real*8 function derf, provided c by clint cummins for use in quail5.1 . the modification using c the fact that phee(z)=0.5+0.5*erf(z/sqr2) was done by vah . c c specifications for arguments c sqr2, pt5 added by vah c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ double precision derf,y,sqr2,pt5 c specifications for local variables dimension p(5),q(4),p1(9),q1(8),p2(6),q2(5) double precision p,q,p1,q1,p2,q2,xmin,xlarge,sqrpi,x, * res,xsq,xnum,xden,xi integer isw,i c added by vah data sqr2 /1.41421356237309d00/,pt5/0.5d0/ c coefficients for 0.0 .le. y .lt. c .477 data p(1)/113.8641541510502d0/, * p(2)/377.4852376853020d0/, * p(3)/3209.377589138469d0/, * p(4)/.1857777061846032d0/, * p(5)/3.161123743870566d0/ data q(1)/244.0246379344442d0/, * q(2)/1282.616526077372d0/, * q(3)/2844.236833439171d0/, * q(4)/23.60129095234412d0/ c coefficients for .477 .le. y c .le. 4.0 data p1(1)/8.883149794388376d0/, * p1(2)/66.11919063714163d0/, * p1(3)/298.6351381974001d0/, * p1(4)/881.9522212417691d0/, * p1(5)/1712.047612634071d0/, * p1(6)/2051.078377826071d0/, * p1(7)/1230.339354797997d0/, * p1(8)/2.153115354744038d-8/, * p1(9)/.5641884969886701d0/ data q1(1)/117.6939508913125d0/, * q1(2)/537.1811018620099d0/, * q1(3)/1621.389574566690d0/, * q1(4)/3290.799235733460d0/, * q1(5)/4362.619090143247d0/, * q1(6)/3439.367674143722d0/, * q1(7)/1230.339354803749d0/, * q1(8)/15.74492611070983d0/ c coefficients for 4.0 .lt. y data p2(1)/-3.603448999498044d-01/, * p2(2)/-1.257817261112292d-01/, * p2(3)/-1.608378514874228d-02/, * p2(4)/-6.587491615298378d-04/, * p2(5)/-1.631538713730210d-02/, * p2(6)/-3.053266349612323d-01/ data q2(1)/1.872952849923460d0/, * q2(2)/5.279051029514284d-01/, * q2(3)/6.051834131244132d-02/, * q2(4)/2.335204976268692d-03/, * q2(5)/2.568520192289822d0/ c constants data xmin/1.0d-10/,xlarge/6.375d0/ data sqrpi/.5641895835477563d0/ c first executable statement c modified by vah x = y/sqr2 isw = 1 if (x.ge.0.0d0) go to 5 isw = -1 x = -x 5 if (x.lt..477d0) go to 10 if (x.le.4.0d0) go to 25 if (x.lt.xlarge) go to 35 res = 1.0d0 go to 50 c abs(y) .lt. .477, evaluate c approximation for derf 10 if (x.lt.xmin) go to 20 xsq = x*x xnum = p(4)*xsq+p(5) xden = xsq+q(4) do 15 i=1,3 xnum = xnum*xsq+p(i) xden = xden*xsq+q(i) 15 continue res = x*xnum/xden go to 50 20 res = x*p(3)/q(3) go to 50 c .477 .le. abs(y) .le. 4.0 c evaluate approximation for derf 25 xsq = x*x xnum = p1(8)*x+p1(9) xden = x+q1(8) do 30 i=1,7 xnum = xnum*x+p1(i) xden = xden*x+q1(i) 30 continue res = xnum/xden go to 45 c 4.0 .lt. abs(y), evaluate c minimax approximation for derf 35 xsq = x*x xi = 1.0d0/xsq xnum = p2(5)*xi+p2(6) xden = xi+q2(5) do 40 i=1,4 xnum = xnum*xi+p2(i) xden = xden*xi+q2(i) 40 continue res = (sqrpi+xi*xnum/xden)/x 45 res = res*dexp(-xsq) res = 1.0d0-res 50 if (isw.eq.-1) res = -res derf = res c added by vah pheev=pt5+pt5*derf return end c c c double precision function pheev2(z1,z2,r) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c program to compute bivariate normal cdf c cacm algorith number 462 c based on algorithm in owen, annals of math stat, 1956:1075-90 c accuracy is controlled by idig--number of sig digits to right of c decimal point in answer. more accurate than old hausman routine c with idig=15 time is approximately .05 seconds on pc w/8087 c time varies with values. can be faster c c modified by vah 6:49 pm, tuesday, december 9, 1986 c c routines used: c pheev(z) is standard cumulative normal at z c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision(a-h,o-z) implicit integer*2(i-n) external pheev ccc data twopi /6.2831853071795865770d00/ ccc data pt5,pt25,zero,one,two /.5d0,.25d0,0.d0,1.d0,2.d0/ if(r .eq. 0.d0)then pheev2=pheev(z1)*pheev(z2) return endif sqr=1.d0-r**2 if(sqr .le. 0.d0)then write(*,801) 801 format(' rho on boundary or outside range: rho=',d13.5) stop endif sqr=dsqrt(sqr) ah=-z1 ak=-z2 twopi=6.283185307179587d0 b=0.d0 ccc idig=15 gh=.5d0*pheev(-ah) gk=.5d0*pheev(-ak) ccc if(idig .ne. 15)then ccc c0n=.5d0*twopi/(10.d0**idig) ccc else ccc c0n=twopi*.5d-15 c0n=twopi*.5d-13 ccc endif if(ah) 170,150,170 150 if(ak) 190,160,190 160 b=datan(r/sqr)/twopi+.25d0 go to 350 170 b=gh if(ah*ak) 180,200,190 180 b=b-.5d0 190 b=b+gk if(ah)200,340,200 200 wh=-ah wk=(ak/ah-r)/sqr gw=2.d0*gh is=-1 210 sgn=-1.d0 t=0.d0 if(wk) 220,320,220 220 if(dabs(wk)-1.d0)270,230,240 230 t=wk*gw*(1.d0-gw)*.5d0 go to 310 240 sgn=-sgn wh=wh*wk g2=pheev(wh) wk=1.d0/wk if(wk)250,260,260 250 b=b+.5d0 260 b=b-(gw+g2)*.5d0+gw*g2 270 h2=wh*wh a2=wk*wk h4=h2*.5d0 ex=dexp(-h4) w2=h4*ex ap=1.d0 s2=ap-ex sp=ap s1=0.d0 sn=s1 c0nex=dabs(c0n/wk) go to 290 280 sn=sp sp=sp+1.d0 s2=s2-w2 w2=w2*h4/sp ap=-ap*a2 290 cn=ap*s2/(sn+sp) s1=s1+cn if(dabs(cn)-c0nex)300,300,280 300 t=(datan(wk)-wk*s1)/twopi 310 b=b+sgn*t 320 if(is)330,350,350 330 if(ak)340,350,340 340 wh=-ak wk=(ah/ak-r)/sqr gw=2.d0*gk is=1 go to 210 350 if(b .lt. 0.d0)then pheev2=0.d0 elseif(b .gt. 1.d0)then pheev2=1.d0 else pheev2=b endif return end subroutine pheev3(z,x1,x2,x3,rho12,rho31,rho23,*) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c vah 12:55 am, sunday, may 22, 1988 c based on farber's program tnorm c c for the trivariate normal r.v. (x,y,w) with zero means, unit c variance and corr(x,y)=rho12,corr(y,w)rho23,corr(w,x)=rho31,tnor c uses the method of section 2 of steck (1958) to compute z= pr c (x0 as required by inverse interpolation: c***** sfeex0=osr2pi*dexp(pt5neg*x0**2) t=(p-pheev(x0))/sfeex0 ph=x0+t+x0*t**2/two+(two*x0**2+one)*t**3/six c***** c back to original x0<0, ph<0: c***** ph=-ph if(iflip.eq.0) ph=-ph phinv=ph return 998 phinv=18.d0 c*****a crude way of dealing with ppp=1 return 999 phinv=-18.d0 c*****a crude way of dealing with ppp=0 return end c c c double precision function randu(iseed) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c vah : 28 mar 86 c 1:07 pm, monday, june 30, 1986 c***** c*****this is a pseudo random number generator, c*****randu being the realization of a uniform on [0,1] c*****pseudo r.v.. requires a seed (iseed) which is updated c*****and returned. written by v.a. hajivassiliou, 24sep84, c*****to exactly reproduce the pr1me uniform pseudo random c*****generator rand$a(iseed). c*****ref: c*****lewis,goodman & miller,ibm systems journal,v.8#2 1969,pp.136-145 c***** c*****uses multiplier number 1 for fishman and moore, "a statistical c*****evaluation of multiplicative congruential random number generators c*****with modulus 2**31-1", jasa 77, march 1982. this was found to c*****be quite acceptable in the tests, and the fastest. c c***** c*****note: the imsl function ggubs is exactly the same algorithm, c***** except that it divides by 2147483648 , potentially c***** speeding the process up on 32 bit machines c***** ccc can't re-mention randu in rm profort: ccc double precision randu,xm,xb,seed,xnews c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit integer*4 (i-n) double precision xm,xb,seed,xnews integer*4 iseed data xb,xm/16807.d0,2147483647.d0/ seed=dble(iseed) xnews=dmod(xb*seed,xm) randu=xnews/xm ccc iseed=intl(xnews) ccc for rmpf: iseed=int(xnews) ccc iseed=jfix(xnews) return end double precision function ranorm(xm,sd,iseed) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c********************************************************************** c***** random number generators package c***** vassilis hajivassiliou october 1984 c***** last modified : 5:45 pm, wednesday, september 24, 1986 c********************************************************************** c this set contains: c ***************** c 1. ranorm(xm,sd,iseed) : (tsp) c 2. binrmb(x,y,xm,ym,xs,ys,rho,iseed) c 3. birnrm(x,y,xm,ym,xs,ys,rho,iseed) c 4. phinv(p) c 5. randb(prob,iseed) c 6. randc(xm,scale,iseed) c 7. randcv(xm,scale,iseed) c 8. rande(xm,iseed) c 9. randl(xm,sd,iseed) c 10. randp(xlambda,iseed) c 11. randpa(xtheta,iseed) c 12. randu(iseed) c 13. ranorv(xm,sd,iseed) : (using phinv) c 14. ranrmb(xm,sd,iseed) c 15. randde(xm,sd,iseed) c c********************************************************************** implicit double precision(a-h,o-z) integer*4 iseed c********************************************************************** c c normal(xm,sd) random number generator . vassilis hajivassiliou c c october 1984 c c this is from cacm algorithm #488, originally implemented c for use in tsp4.0 by bronwyn h. hall(oct82). c c********************************************************************** c subprogram(s) used: randu c********************************************************************** c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ dimension d(60),d1(30),d2(30) equivalence(d1(1),d(1)),(d2(1),d(31)) c data d1/ * 0.674489750d0,0.475859630d0,0.383771164d0,0.328611323d0, * 0.291142827d0,0.263684322d0,0.242508452d0,0.225567444d0, * 0.211634166d0,0.199924267d0,0.189910758d0,0.181225181d0, * 0.173601400d0,0.166841909d0,0.160796729d0,0.155349717d0, * 0.150409384d0,0.145902577d0,0.141770033d0,0.137963174d0, * 0.134441762d0,0.131172150d0,0.128125965d0,0.125279090d0, * 0.122610883d0,0.120103560d0,0.117741707d0,0.115511892d0, * 0.113402349d0,0.111402720d0/ data d2/ * 0.109503852d0,0.107697617d0,0.105976772d0,0.104334841d0, * 0.102766012d0,0.101265052d0,0.099827234d0,0.098448282d0, * 0.097124309d0,0.095851778d0,0.094627461d0,0.093448407d0, * 0.092311909d0,0.091215482d0,0.090156838d0,0.089133867d0, * 0.088144619d0,0.087187293d0,0.086260215d0,0.085361834d0, * 0.084490706d0,0.083645487d0,0.082824924d0,0.082027847d0, * 0.081253162d0,0.080499844d0,0.079766932d0,0.079053527d0, * 0.078358781d0,0.077681899d0/ c u = randu(iseed) a = 0.0d0 i = 0 10 continue if (u.eq.1.0d0) u = randu(iseed) u = u+u if (u.lt.1.0d0) go to 20 u = u-1.0d0 i = i+1 a = a-d(i) go to 10 20 continue w = d(i+1)*u v = w*(0.5d0*w-a) 30 continue u = randu(iseed) if (v.le.u) go to 40 v = randu(iseed) if (u.gt.v) go to 30 u = (v-u)/(1.0d0-u) go to 20 40 continue u = (u-v)/(1.0d0-v) u = u+u if (u.lt.1.0d0) go to 50 u = u-1.0d0 rn = w-a go to 100 50 continue rn = a-w 100 continue ranorm=xm+sd*rn return end double precision function sfee(x) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision(a-h,o-z) data osr2pi/0.39894228040143267793d00/,pt5neg/-0.5d0/ sfee=0.d0 xsq=x*x if(xsq.gt.1000.d0) return sfee=osr2pi*dexp(pt5neg*xsq) return end subroutine sorts(y,n) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c ccc subroutine sorts(y,n,ierr) c ccc integer*4 n,ierr integer*4 n double precision y(n) c c shell sort n values in t() from smallest to largest. c ***** c c note that local system sort utilities are likely to be c more efficient, and should be substituted whenever possible. c c "abc's of eda" manual 1981 p.313 c written 11/03/85 17:17:29 c c local variables c ccc integer*2 i, j, j1, igap, nmg c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ integer*4 i, j, j1, igap, nmg double precision temp c if(n .ge. 1) goto 10 ccc ierr=1 goto 999 10 if(n .eq.1) goto 999 c c one element is always sorted c igap=n 20 igap=igap/2 nmg=n-igap do 40 j1=1,nmg i=j1+igap c c do j=j1,1,-igap c j=j1 30 if(y(j) .le. y(i)) goto 40 c c swap out-of-order pair c temp=y(i) y(i)=y(j) y(j)=temp c c keep old pointer for next time through c i=j j=j-igap if(j .ge. 1) goto 30 40 continue if(igap .gt. 1) goto 20 999 return end c c c subroutine srtrnk(yvec,n,ivec) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c vah 8:11 pm, tuesday, october 28, 1986 c 16:50:51.97, sun, apr 9, 1989 c implicit double precision(a-h,o-z) implicit integer (i-n) c c shell sort n values in yvec() from smallest to largest. c ***** c vah 7:29 pm, tuesday, october 28, 1986 c modified to return also array ivec. c on entry, ivec(1) controls the operations to be performed: c if ivec(1) ge 0: yvec is sorted from min to max (call this xsort) c ivec will contain the original locations of the c elements of yvec (call this iloc) c if ivec(1) lt 0: yvec is sorted from min to max in place, i.e. c the original yvec is returned c ivec will contain the ranks of the yvec c (call this irnk) c----------------------------------------------------------------------- c the full relations are: c c information returned with information returned with c original ivec(1) ge 0 with ivec(1) lt 0 c -------- ------------------------- ------------------------- c yvec yvec ivec yvec ivec c ---- ---- ---- ---- ---- c x xsort iloc x irnk c - ----- ---- - ---- c 2 1 4 2 2 c 4 2 1 4 3 c 7 4 2 7 5 c 1 5 5 1 1 c 5 7 3 5 4 c c to obtain x and irnk to obtain xsort and iloc c from xsort, iloc from x, irnk c -------------------- ------------------------ c x(iloc(j)=xsort(j) xsort(irnk(j))=x(j) c irnk(iloc(j))=j iloc(irnk(j))=j c c----------------------------------------------------------------------- c c note that local system sort utilities are likely to be c more efficient, and should be substituted whenever possible. c c "abc's of eda" manual 1981 p.313 c written 11/03/85 17:17:29 c c local variables c ccc integer*2 i, j, j1, igap, nmg c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ dimension yvec(n),ivec(n) if(ivec(1).lt.0) then isort=-1 else isort=1 endif c if(n .ge. 1) goto 10 ccc ierr=-1 ivec(1)=1 goto 999 10 if(n .eq.1) goto 999 c c one element is always sorted c do 555 i=1,n 555 ivec(i)=i igap=n 20 igap=igap/2 nmg=n-igap do 40 j1=1,nmg i=j1+igap c c do j=j1,1,-igap c j=j1 30 if(yvec(j) .le. yvec(i)) goto 40 c c swap out-of-order pair c temp=yvec(i) itemp=ivec(i) yvec(i)=yvec(j) ivec(i)=ivec(j) yvec(j)=temp ivec(j)=itemp c c keep old pointer for next time through c i=j j=j-igap if(j .ge. 1) goto 30 40 continue if(igap .gt. 1) goto 20 ccc do 777 j=1,n ccc irnk(ivec(j))=j ccc 777 continue if(isort.le.0) then next0=1 inext=ivec(1) temp0=yvec(1) 887 continue do 888 jjj=1,n c print *,'switching',inext temp1=yvec(inext) yvec(inext)=temp0 temp0=temp1 next=ivec(inext) ivec(inext)=-next0 next0=inext inext=next if(inext.lt.0) then c print *,'inext lt 0 next0',inext,next0 do 886 iii=1,n if(ivec(iii).gt.0) then inext=ivec(iii) next0=iii c print *,'restarting with inext,next0:',inext,next0 goto 887 endif 886 continue goto 889 endif 888 continue 889 continue do 884 iii=1,n ivec(iii)=-ivec(iii) 884 continue endif c print *,'ivec',(ivec(kkk),kkk=1,n) c print *,'yvec',(yvec(kkk),kkk=1,n) 999 return end subroutine uc(ch) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c changes the character variable ch to upper case c vah 5:13 pm, thursday, april 3, 1986 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ integer*4 length character*(*) ch character*1 c1 length=len(ch) do 1 i=1,length c1=ch(i:i) ic=ichar(c1) c print *,'ic:',ic if(ic.gt.96.and.ic.lt.123) ch(i:i)=char(ic-32) 1 continue return end subroutine when(nfunit) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ integer*4 nfunit integer*4 itim(5),idat(3) character*3 cm(12) cm(1)='jan' cm(2)='feb' cm(3)='mar' cm(4)='apr' cm(5)='may' cm(6)='jun' cm(7)='jul' cm(8)='aug' cm(9)='sep' cm(10)='oct' cm(11)='nov' cm(12)='dec' call timdat(itim,idat) if(idat(3).lt.100) idat(3)=idat(3)+1900 write(nfunit,100) idat(1),cm(idat(2)),idat(3),(itim(i),i=1,4) 100 format('0','date:',i4,' ',a3,' ',i4,2x, * 'time:',2x,i2.2,':',i2.2,':',i2.2,'.',i2.2/ * 1x,'----', 16x,'----') return end