V3FIT
getspline.f
1  SUBROUTINE getspline(amat, splnots, hk, delse, hs, indexs, isort,
2  1 ndata0, nots)
3  USE vsvd0
4  USE vparams, ONLY: zero, epstan
5  IMPLICIT NONE
6 C-----------------------------------------------
7 C D u m m y A r g u m e n t s
8 C-----------------------------------------------
9  INTEGER ndata0, nots
10  REAL(rprec) hs
11  INTEGER, DIMENSION(ndata0) :: indexs, isort
12  REAL(rprec), DIMENSION(nots,ndata0) :: amat
13  REAL(rprec), DIMENSION(nots) :: splnots, hk
14  REAL(rprec), DIMENSION(ndata0) :: delse
15 C-----------------------------------------------
16 C L o c a l V a r i a b l e s
17 C-----------------------------------------------
18  INTEGER, DIMENSION(nots) :: nk
19  INTEGER :: i, js, ia, k, k1, ib, j, nb
20  REAL(rprec), DIMENSION(ndata0) :: w, w1, u, u1, snodes
21  REAL(rprec), DIMENSION(nots,ndata0) :: bmat
22 C-----------------------------------------------
23 
24 !
25 ! ON EXIT, AMAT = AMAT + BMAT, WHERE BMAT IS THE 2nd
26 ! DERIVATIVE COEFFICIENT MATRIX ARRAY, MULTIPLIED BY JACOBIAN
27 ! AND AMAT (ON RHS) WAS FUNCTION COEFFICIENT MATRIX ARRAY
28 !
29 
30 !
31 ! SORT KNOT POSITIONS IN ASCENDING ORDER IN S-SPACE
32 ! USE SQRT(S) KNOT POSITIONS FOR IMPROVED RESOLUTION
33 ! NOTE: SNODES(I) IS THE VALUE OF SQRT(S) CORRESPONDING TO
34 ! THE MESH VALUES CORRESPONDING TO DELSE, INDEXS (COMPUTED OUTSIDE
35 ! THIS PROGRAM)
36 !
37 
38  DO i = 1, ndata0
39  js = indexs(i)
40  snodes(i) = sqrt(hs*((js - 1) + delse(i)))
41 ! IF (snodes(i) .le. zero) snodes(i) = epstan
42  END DO
43 
44 ! Avoid roundoff error in SPLININT
45  IF( snodes(ndata0) .gt. splnots(nots) )
46  1 snodes(ndata0) = splnots(nots)
47 
48  CALL sort_data (snodes, isort, ndata0)
49 
50 !
51 ! COMPUTE MATRIX COEFFICIENTS RELATING SPLINE AT SPLNOTS
52 ! TO REAL-SPACE FUNCTION AT SORTED MESH POINTS RMESH
53 !
54  amat(:nots,:ndata0) = zero
55  bmat(:nots,:ndata0) = zero
56 
57 !
58 ! SETUP SPLINE PARAMETERS AT EACH TIME STEP, SINCE SNODES
59 ! MAY BE CHANGING DYNAMICALLY IN S-SPACE
60 !
61  CALL setup_int(splnots,snodes,hk,w,w1,u,u1,nk,nots,ndata0)
62 
63  ia = 1
64  DO k = 1, nots - 1
65  IF (nk(k) .gt. 0) THEN
66  k1 = k + 1
67  ib = ia + nk(k) - 1
68  amat(k,ia:ib) = amat(k,ia:ib) + w(ia:ib)
69  bmat(k,ia:ib) = bmat(k,ia:ib) + u(ia:ib)
70  amat(k1,ia:ib) = amat(k1,ia:ib) + w1(ia:ib)
71  bmat(k1,ia:ib) = bmat(k1,ia:ib) + u1(ia:ib)
72  ia = ib + 1
73  ENDIF
74  END DO
75 
76  IF (ib .ne. ndata0) stop 'ib!=ndat'
77  nb = ideriv
78 
79  DO j = 1, ndata0
80  bmat(nots,j) = 0.
81  CALL jacprod (bmat(1,j), hk, nots, nb)
82  amat(:nots,j) = amat(:nots,j) + bmat(:nots,j)
83  END DO
84 
85  END SUBROUTINE getspline