Click tab to download
! BSPLINES ! Purpose: evaluate B-spline basis functions ! linear, quadratic or cubic ! equidistant (default) or at user-specified knots ! Input: via command line options (-h will list them) or interactively ! program will prompt for esential options not given on command line; ! ages and knot positions can be given in separate files ! Choose ! degree: command line option -d1, -d2 or -d3 for linear, quadratic or cubic; ! if not given, program will prompt for it ! "ages": command line option -a to read min. & max. age and step size ! interactively; if not given program will prompt for it, ! if given, program will expect to read ages from file "agesfile", ! default name "ages.dat" ! knots: command line option -eqx specificies a spline function with equal ! intervals, defined by x knots (x must be integer); if not given, ! program will prompt for the no. of knots. ! to select a spline with user-selected knots, a "knotsfile" (default ! name "knots" must exist in the working directory; program will ! count no. of knots in the file ! Format of "agesfile": ! one age per line ! each line comprised of two space-separated values: running no. ! (integer) and value for age (integer) ! Format for "knotsfile": ! one line per knot containing a single (real) value, ! knots should be in ascending order; ! it is the user's responsibility to ensure that knots are sensible & ! that the external knots & the range of ages specified are consistent! ! Output file "bspfile" with default name "bspline_Xn.dat" where X is the ! degree of the spline (L, Q or C) and n is the no. of knots given: ! lines in order of ages ! one line per age with n coefficients (real), followed by value for ! age; adjust format for output as necessary ! Compilation note: 'gfortran' may not accommadate the intrinsic functions ! 'iargc' and 'getarg' in subroutine cmdline; replace by ! 'k = command_argument_count()' and ! 'call get_command_argument( iarg, arg )' respectively ! ====================================================================== program Bsplines ! ===================================================@C-km-2004-======== implicit none ! -------------------- change as desired ---------------------------------- ! default name for file with user-chosen knots character(len=20) :: knotsfile ='knots ' ! default name for file with "ages" at which to evaluate B-spline function character(len=20) :: agesfile = 'ages.dat ' ! default name for output file character(len=16) :: bsplfile = 'bspline_X00.dat' ! format for output file character(len=12) :: fmt = '(00f12.8,i8)' ! ------------------- end of changes -------------------------------------- character(len=1), dimension(3) :: degree = (/ 'L', 'Q', 'C' /) real(8), dimension(:), allocatable :: kvec, tvec, w, p integer, dimension(:), allocatable :: iiage real(8), dimension(:,:,:), allocatable :: mm integer :: nknot, ndeg, nord, mknot, nage, nrr integer :: ii, i, j, iage, ia, iage1, iage2, & & istep real(8) :: xx, age, age1, age2, aa logical :: verbose, degset, kntset, equal, & & ages, lex ! read command line options call cmdline ! degree of basis functions if(.not.degset)then write(*,*)'degree ? (1=linear, 2=quad., 3=cubic)' read(*,*) ndeg; if(ndeg <1 .or. ndeg >3 ) stop 'wrong degree' end if nord = ndeg+1 ! order of fit ! ages for which to evaluate basis functions 100 continue if( ages ) then write(*,*)'Minimum age ?'; read(*,*)iage1 write(*,*)'Maximum age ?'; read(*,*)iage2 if(iage2.le.iage1)then write(*,*)'error : max. age must be greater than min. age' stop 'Bsplines' end if write(*,*)'Step size ?'; read(*,*) istep allocate( iiage(iage2-iage1+1) ); nage=0 do iage=iage1,iage2,istep nage=nage+1; iiage(nage) = iage end do if(iiage(nage)<iage2)write(*,*) & & 'Warning : Age range not a multiple of step size !' else inquire(file = agesfile,exist=lex) if( .not. lex)then write(*,*)'File',agesfile,'does not exist !' ages=.true.; go to 100 end if open(1,file=agesfile,status='old', action='read' ) nage = 0 do read(1,*,iostat=ii) xx; if(ii /= 0) exit nage = nage + 1 end do allocate( iiage(nage)); rewind(1) do read(1,*,iostat=ii) nage, iiage(nage); if(ii/=0)exit if(verbose) write(*,*)'age', nage, iiage(nage) end do close(1) end if if( verbose ) then write(*,'(a,t30,1h=,i8)')'No. of ages',nage write(*,'(a,t30,1h=,2i8)')'Age range',iiage(1),iiage(nage) end if ! where are the knots ? 200 continue if( equal ) then if( .not. kntset ) then write(*,'(a)')'No. of knots ?'; read(*,*) nknot; kntset=.true. end if age1=iiage(1); age2=iiage(nage); aa= ( age2 - age1 ) / dble( nknot - 1 ) allocate( kvec(nknot) ) do i = 1, nknot kvec(i) = age1 + (i-1)*aa end do mknot = nknot + 2*ndeg; nrr = nknot + ndeg -1 allocate( tvec(mknot) ) tvec(nord:nord+nknot-1) = kvec do i=1,ndeg tvec(i)= age1 - (ndeg-i+1)*aa end do do i=1,ndeg tvec(nknot+nord-1+i) = age2 + i* aa end do else inquire(file=knotsfile,exist=lex) if(.not.lex) then write(*,'(a,1h",a,1h",a)')'File ',trim(knotsfile),' does not exist !' equal=.true.; go to 200 end if open(1,file=knotsfile,status='old',action='read') ; nknot=0 do read(1,*,iostat=ii)xx; if(ii /= 0)exit nknot = nknot+1 end do write(*,*) 'no. of user-specified knots =', nknot allocate( kvec(nknot) ); rewind(1) do i=1,nknot read(1,*,iostat=ii)kvec(i); if(ii /= 0)exit if(verbose)write(*,*)'knot',i,kvec(i) end do close(1) if(iiage(1)<kvec(1) .or. iiage(nage) .ge. kvec(nknot) ) then ! check that range of ages & external knots are consistent write(*,*)'range of ages and knots specified disagree !' write(*,*)'external knots given =',kvec(1),kvec(nknot) write(*,*)'range of ages found =',iiage(1),iiage(nage) stop 'Bsplines' end if mknot = nknot + 2*ndeg; nrr = nknot + ndeg -1 ! no. of RR coefficients allocate( tvec(mknot) ) tvec(1:ndeg)=kvec(1); tvec(nord:nord+nknot-1) = kvec tvec(nord+nknot:mknot) = kvec(nknot) end if if(verbose)write(*,*)'(Extended) Vector of knots ',tvec ! basis functions of order 1 allocate( mm(mknot,0:ndeg,nage) ) mm = 0.d0 do i = 1,nage iage = iiage(i) do j = 1,mknot-1 if(tvec(j) == tvec(j+1))cycle if(iage.ge.tvec(j).and.iage.lt.tvec(j+1)) mm(j,0,i)=1.d0 end do end do ! higher order basis functions do ii = 1, ndeg do ia = 1,nage age = iiage(ia) do j = 1,mknot-ii if( tvec(j+ii) > tvec(j) ) mm(j,ii,ia) = & & ( age - tvec(j) ) / ( tvec(j+ii)-tvec(j)) * mm(j,ii-1,ia) if(j+ii+1 > mknot) cycle if( tvec(j+ii+1) > tvec(j+1) ) mm(j,ii,ia) = mm(j,ii,ia) + & & ( tvec(j+ii+1) - age) / ( tvec(j+ii+1) - tvec(j+1)) * mm(j+1,ii-1,ia) end do end do end do ! write out file of covariables bsplfile(9:9)=degree(ndeg) if(nrr <10 )then write(bsplfile(11:11),'(i1)') nrr else write(bsplfile(10:11),'(i2)') nrr end if open(11,file=bsplfile,status='unknown',form='formatted') write( fmt(2:3),'(i2)' ) nrr do ia = 1, nage write(11,fmt) mm(:nrr,ndeg,ia), iiage(ia) end do close(11) ! write out a summary write(*,*)'number of ages = ',nage write(*,*)' ... range ',iiage(1),iiage(nage) write(*,*)'degree of polynomial = ',ndeg write(*,*)'number of knots = ',nknot write(*,*)'number of RR coefficients = ',nrr if( equal ) write(*,*)'Uniform spline with auxiliary knots outside range' write(*,'(1x,a,1h",a,1h")')'output file : ',bsplfile contains ! ================== subroutine cmdline ! ================== ! routine to evaluate command line options integer :: iargc, i, k, l character(len=30) :: arg external getarg ! default value for options verbose =.false.; degset=.false.; kntset=.false.; equal=.false.; ages =.false. k = iargc() ! no. of arguments do i = 1, k call getarg (i,arg); arg = adjustl(arg) if(arg(1:2) == '-v' .or. arg(1:2) == '-V') then verbose = .true. else if(arg(1:3) == '-eq' .or. arg(1:3) == '-EQ' ) then equal=.true.; read(arg(4:),*)nknot; kntset=.true. else if(arg(1:2) == '-a' .or. arg(1:2) == '-A')then ages = .true. else if( arg(1:2) == '-d' .or. arg(1:2) == '-D') then read(arg(3:),*) ndeg; if(ndeg <1 .or. ndeg >3 ) stop 'wrong degree' degset = .true. else if( arg(1:2) =='-h' .or. arg(1:2) == '-H') then write(*,*)'Valid command line options :' write(*,'(a,t10,a)')'-h','produces this message' write(*,'(a,t10,a)')'-u','how to use it' write(*,'(a,t10,a)')'-v','verbose screen output' write(*,'(a,t10,a)')'-a','specify range of ages & step size' write(*,'(a,t10,a)')'-dx','for x=1,2,3 : degree of spline' write(*,'(a,t10,a)')'-eqx','for x=1,2,... : uniform spline with x knots' else if( arg(1:2) == '-u' .or. arg(1:2) == '-U' ) then write(*,*)'Program BSPLINES : set up B-spline functions ' write(*,*)'------------------------------------------------' write(*,*)'Input: file "ages.out"' write(*,*)'Input: File "knots" with user-defined knots (one per line)' write(*,*)'Input: Options via command line or interactively' write(*,*)' - degree of spline (1=linear, 2=quadratic, 3=cubic)' write(*,*)' - absence of "knots" selects a uniform spline with' write(*,*)' auxiliary knots equi-distant outside range' write(*,*)' - for "ages.out" absent, user is prompted for age' write(*,*)' range and step size '; stop end if end do if( verbose ) then write(*,*)' **** Command line options ***' if(degset) write(*,'(1x,a,t30,1h=,i12)')'Degree of spline',ndeg if(equal) write(*,'(1x,a,t30,1h=,i12,1x,a)')'Uniform spline with', & & nknot,'knots' end if return end subroutine cmdline end program Bsplines