[1765] | 1 | MODULE poly_m |
---|
| 2 | |
---|
| 3 | IMPLICIT NONE |
---|
| 4 | |
---|
| 5 | INTEGER, PARAMETER, private :: NPAR_POLY=8 |
---|
| 6 | |
---|
| 7 | INTERFACE poly |
---|
| 8 | MODULE PROCEDURE poly_rr,poly_rrv,poly_dd,poly_ddv,& |
---|
| 9 | poly_rc,poly_cc,poly_msk_rrv,poly_msk_ddv |
---|
| 10 | END INTERFACE |
---|
| 11 | |
---|
| 12 | private poly_rr,poly_rrv,poly_dd,poly_ddv,& |
---|
| 13 | poly_rc,poly_cc,poly_msk_rrv,poly_msk_ddv |
---|
| 14 | |
---|
| 15 | CONTAINS |
---|
| 16 | |
---|
| 17 | FUNCTION poly_rr(x,coeffs) |
---|
| 18 | REAL, INTENT(IN) :: x |
---|
| 19 | REAL, DIMENSION(:), INTENT(IN) :: coeffs |
---|
| 20 | REAL :: poly_rr |
---|
| 21 | REAL :: pow |
---|
| 22 | REAL, DIMENSION(:), ALLOCATABLE :: vec |
---|
| 23 | INTEGER :: i,n,nn |
---|
| 24 | n=size(coeffs) |
---|
| 25 | if (n <= 0) then |
---|
| 26 | poly_rr=0.0 |
---|
| 27 | else if (n < NPAR_POLY) then |
---|
| 28 | poly_rr=coeffs(n) |
---|
| 29 | do i=n-1,1,-1 |
---|
| 30 | poly_rr=x*poly_rr+coeffs(i) |
---|
| 31 | end do |
---|
| 32 | else |
---|
| 33 | allocate(vec(n+1)) |
---|
| 34 | pow=x |
---|
| 35 | vec(1:n)=coeffs |
---|
| 36 | do |
---|
| 37 | vec(n+1)=0.0 |
---|
| 38 | nn=ishft(n+1,-1) |
---|
| 39 | vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2) |
---|
| 40 | if (nn == 1) exit |
---|
| 41 | pow=pow*pow |
---|
| 42 | n=nn |
---|
| 43 | end do |
---|
| 44 | poly_rr=vec(1) |
---|
| 45 | deallocate(vec) |
---|
| 46 | end if |
---|
| 47 | END FUNCTION poly_rr |
---|
| 48 | !BL |
---|
| 49 | FUNCTION poly_dd(x,coeffs) |
---|
| 50 | DOUBLE PRECISION, INTENT(IN) :: x |
---|
| 51 | DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: coeffs |
---|
| 52 | DOUBLE PRECISION :: poly_dd |
---|
| 53 | DOUBLE PRECISION :: pow |
---|
| 54 | DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: vec |
---|
| 55 | INTEGER :: i,n,nn |
---|
| 56 | n=size(coeffs) |
---|
| 57 | if (n <= 0) then |
---|
| 58 | poly_dd=0.0d0 |
---|
| 59 | else if (n < NPAR_POLY) then |
---|
| 60 | poly_dd=coeffs(n) |
---|
| 61 | do i=n-1,1,-1 |
---|
| 62 | poly_dd=x*poly_dd+coeffs(i) |
---|
| 63 | end do |
---|
| 64 | else |
---|
| 65 | allocate(vec(n+1)) |
---|
| 66 | pow=x |
---|
| 67 | vec(1:n)=coeffs |
---|
| 68 | do |
---|
| 69 | vec(n+1)=0.0d0 |
---|
| 70 | nn=ishft(n+1,-1) |
---|
| 71 | vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2) |
---|
| 72 | if (nn == 1) exit |
---|
| 73 | pow=pow*pow |
---|
| 74 | n=nn |
---|
| 75 | end do |
---|
| 76 | poly_dd=vec(1) |
---|
| 77 | deallocate(vec) |
---|
| 78 | end if |
---|
| 79 | END FUNCTION poly_dd |
---|
| 80 | !BL |
---|
| 81 | FUNCTION poly_rc(x,coeffs) |
---|
| 82 | COMPLEX, INTENT(IN) :: x |
---|
| 83 | REAL, DIMENSION(:), INTENT(IN) :: coeffs |
---|
| 84 | COMPLEX :: poly_rc |
---|
| 85 | COMPLEX :: pow |
---|
| 86 | COMPLEX, DIMENSION(:), ALLOCATABLE :: vec |
---|
| 87 | INTEGER :: i,n,nn |
---|
| 88 | n=size(coeffs) |
---|
| 89 | if (n <= 0) then |
---|
| 90 | poly_rc=0.0 |
---|
| 91 | else if (n < NPAR_POLY) then |
---|
| 92 | poly_rc=coeffs(n) |
---|
| 93 | do i=n-1,1,-1 |
---|
| 94 | poly_rc=x*poly_rc+coeffs(i) |
---|
| 95 | end do |
---|
| 96 | else |
---|
| 97 | allocate(vec(n+1)) |
---|
| 98 | pow=x |
---|
| 99 | vec(1:n)=coeffs |
---|
| 100 | do |
---|
| 101 | vec(n+1)=0.0 |
---|
| 102 | nn=ishft(n+1,-1) |
---|
| 103 | vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2) |
---|
| 104 | if (nn == 1) exit |
---|
| 105 | pow=pow*pow |
---|
| 106 | n=nn |
---|
| 107 | end do |
---|
| 108 | poly_rc=vec(1) |
---|
| 109 | deallocate(vec) |
---|
| 110 | end if |
---|
| 111 | END FUNCTION poly_rc |
---|
| 112 | !BL |
---|
| 113 | FUNCTION poly_cc(x,coeffs) |
---|
| 114 | COMPLEX, INTENT(IN) :: x |
---|
| 115 | COMPLEX, DIMENSION(:), INTENT(IN) :: coeffs |
---|
| 116 | COMPLEX :: poly_cc |
---|
| 117 | COMPLEX :: pow |
---|
| 118 | COMPLEX, DIMENSION(:), ALLOCATABLE :: vec |
---|
| 119 | INTEGER :: i,n,nn |
---|
| 120 | n=size(coeffs) |
---|
| 121 | if (n <= 0) then |
---|
| 122 | poly_cc=0.0 |
---|
| 123 | else if (n < NPAR_POLY) then |
---|
| 124 | poly_cc=coeffs(n) |
---|
| 125 | do i=n-1,1,-1 |
---|
| 126 | poly_cc=x*poly_cc+coeffs(i) |
---|
| 127 | end do |
---|
| 128 | else |
---|
| 129 | allocate(vec(n+1)) |
---|
| 130 | pow=x |
---|
| 131 | vec(1:n)=coeffs |
---|
| 132 | do |
---|
| 133 | vec(n+1)=0.0 |
---|
| 134 | nn=ishft(n+1,-1) |
---|
| 135 | vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2) |
---|
| 136 | if (nn == 1) exit |
---|
| 137 | pow=pow*pow |
---|
| 138 | n=nn |
---|
| 139 | end do |
---|
| 140 | poly_cc=vec(1) |
---|
| 141 | deallocate(vec) |
---|
| 142 | end if |
---|
| 143 | END FUNCTION poly_cc |
---|
| 144 | !BL |
---|
| 145 | FUNCTION poly_rrv(x,coeffs) |
---|
| 146 | REAL, DIMENSION(:), INTENT(IN) :: coeffs,x |
---|
| 147 | REAL, DIMENSION(size(x)) :: poly_rrv |
---|
| 148 | INTEGER :: i,n,m |
---|
| 149 | m=size(coeffs) |
---|
| 150 | n=size(x) |
---|
| 151 | if (m <= 0) then |
---|
| 152 | poly_rrv=0.0 |
---|
| 153 | else if (m < n .or. m < NPAR_POLY) then |
---|
| 154 | poly_rrv=coeffs(m) |
---|
| 155 | do i=m-1,1,-1 |
---|
| 156 | poly_rrv=x*poly_rrv+coeffs(i) |
---|
| 157 | end do |
---|
| 158 | else |
---|
| 159 | do i=1,n |
---|
| 160 | poly_rrv(i)=poly_rr(x(i),coeffs) |
---|
| 161 | end do |
---|
| 162 | end if |
---|
| 163 | END FUNCTION poly_rrv |
---|
| 164 | !BL |
---|
| 165 | FUNCTION poly_ddv(x,coeffs) |
---|
| 166 | DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: coeffs,x |
---|
| 167 | DOUBLE PRECISION, DIMENSION(size(x)) :: poly_ddv |
---|
| 168 | INTEGER :: i,n,m |
---|
| 169 | m=size(coeffs) |
---|
| 170 | n=size(x) |
---|
| 171 | if (m <= 0) then |
---|
| 172 | poly_ddv=0.0d0 |
---|
| 173 | else if (m < n .or. m < NPAR_POLY) then |
---|
| 174 | poly_ddv=coeffs(m) |
---|
| 175 | do i=m-1,1,-1 |
---|
| 176 | poly_ddv=x*poly_ddv+coeffs(i) |
---|
| 177 | end do |
---|
| 178 | else |
---|
| 179 | do i=1,n |
---|
| 180 | poly_ddv(i)=poly_dd(x(i),coeffs) |
---|
| 181 | end do |
---|
| 182 | end if |
---|
| 183 | END FUNCTION poly_ddv |
---|
| 184 | !BL |
---|
| 185 | FUNCTION poly_msk_rrv(x,coeffs,mask) |
---|
| 186 | REAL, DIMENSION(:), INTENT(IN) :: coeffs,x |
---|
| 187 | LOGICAL, DIMENSION(:), INTENT(IN) :: mask |
---|
| 188 | REAL, DIMENSION(size(x)) :: poly_msk_rrv |
---|
| 189 | poly_msk_rrv=unpack(poly_rrv(pack(x,mask),coeffs),mask,0.0) |
---|
| 190 | END FUNCTION poly_msk_rrv |
---|
| 191 | !BL |
---|
| 192 | FUNCTION poly_msk_ddv(x,coeffs,mask) |
---|
| 193 | DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: coeffs,x |
---|
| 194 | LOGICAL, DIMENSION(:), INTENT(IN) :: mask |
---|
| 195 | DOUBLE PRECISION, DIMENSION(size(x)) :: poly_msk_ddv |
---|
| 196 | poly_msk_ddv=unpack(poly_ddv(pack(x,mask),coeffs),mask,0.0d0) |
---|
| 197 | END FUNCTION poly_msk_ddv |
---|
| 198 | |
---|
| 199 | END MODULE poly_m |
---|