subroutine hpsort_indx(n,da,indx) c c*** Heapsort adapted from Numerical Recipes. da is double precision c*** array which is sorted into ascending order. Indx is an c*** integer array that defines the corresponding permutation and its c*** inverse. That is, c*** c*** da_sort(i) = da_orig(indx(i)) c*** da_orig(indx(i)) = da_sort(i) c*** c*** Analogous relations can be applied to associated arrays that c*** are to be sorted with da as the key. For example, for array db, c*** c*** db_sort(i) = db_orig(indx(i)) c*** db_orig(indx(i)) = db_sort(i) c*** c*** G. Rybicki 3 Feb 95. c INTEGER n DOUBLE PRECISION da(n),datemp INTEGER i,ir,j,k,indx(n),indxtemp c if(n.lt.2)then return endif c do i=1,n indx(i)=i enddo c k=n/2+1 ir=n 10 continue if(k.gt.1)then k=k-1 datemp=da(k) indxtemp=indx(k) else datemp=da(ir) indxtemp=indx(ir) da(ir)=da(1) indx(ir)=indx(1) ir=ir-1 if(ir.eq.1)then da(1)=datemp indx(1)=indxtemp go to 30 endif endif i=k j=k+k 20 if(j.le.ir)then if(j.lt.ir)then if(da(j).lt.da(j+1))then j=j+1 endif endif if(datemp.lt.da(j))then da(i)=da(j) indx(i)=indx(j) i=j j=j+j else j=ir+1 endif go to 20 endif da(i)=datemp indx(i)=indxtemp go to 10 30 continue return end c c***************************************************************************** c subroutine trimult(a,b,c,x,y,n) c c Does the multiplication of a vector by a tridiagonal matrix, i.e., c c y(i) = a(i)*x(i-1) + b(i)*x(i) + c(i)*x(i+1), i=1,...,n, c c with special cases i=1 and i=n, for which the terms containing c a(1) and c(n) are assumed absent. G. Rybicki c INTEGER i,n,n1 DOUBLE PRECISION a(*),b(*),c(*),x(*),y(*) c n1=n-1 y(1)=b(1)*x(1)+c(1)*x(2) do i=2,n1 y(i)=a(i)*x(i-1)+b(i)*x(i)+c(i)*x(i+1) enddo y(n)=a(n)*x(n1)+b(n)*x(n) return end c c***************************************************************************** c subroutine tridiag(a,b,c,x,y,n) c c*** Solves the tridiagonal system c*** c*** a(i)*x(i-1)+b(i)*x(i)+c(i)*x(i+1)=y(i), i=1,n c*** c*** where the terms containing a(1) and c(n) do not appear. c*** Arrays d and e are scratch arrays of dimension NDIM, which c*** must be at least equal to n-1. G. Rybicki c INTEGER i,k,n,n1,NDIM PARAMETER (NDIM=10000) !alter to suit problem double precision a(*),b(*),bb,c(*),d(NDIM),e(NDIM),x(*),y(*) n1=n-1 if(n1.gt.NDIM)then stop 'tridiag: NDIM too small' endif d(1)=-c(1)/b(1) e(1)=y(1)/b(1) do i=2,n1 bb=b(i)+a(i)*d(i-1) d(i)=-c(i)/bb e(i)=(y(i)-a(i)*e(i-1))/bb enddo x(n)=(y(n)-a(n)*e(n1))/(b(n)+a(i)*d(n1)) do k=n1,1,-1 x(k)=e(k)+d(k)*x(k+1) enddo return end c c***************************************************************************** c subroutine tridet(a,b,c,n,logdet,signdet) c c*** Calculate the determinant of the tridiagonal matrix with diagonals c*** a, b, and c. The sign is provided by signdet and the magnitude is c*** provided by its logarithm logdet. The elements a(1) and c(n) c*** are irrelevant and are not used. G. Rybicki c INTEGER i,n,n1 DOUBLE PRECISION a(*),b(*),c(*),logdet,signdet DOUBLE PRECISION bb,dd,ONE PARAMETER (ONE=1.d0) c n1=n-1 signdet=sign(ONE,b(1)) logdet=log(abs(b(1))) dd=-c(1)/b(1) do i=2,n1 bb=b(i)+a(i)*dd dd=-c(i)/bb signdet=signdet*sign(ONE,bb) logdet=logdet+log(abs(bb)) enddo bb=b(n)+a(n)*dd signdet=signdet*sign(ONE,bb) logdet=logdet+log(abs(bb)) return end