![]() |
Millepede-II
V04-00-00
|
00001 00002 ! Code converted using TO_F90 by Alan Miller 00003 ! Date: 2012-03-16 Time: 11:07:10 00004 00022 00034 00035 SUBROUTINE daxpy(n,da,dx,incx,dy,incy) 00036 00037 DOUBLE PRECISION :: dx(*),dy(*),da 00038 INTEGER :: i,incx,incy,ix,iy,m,mp1,n 00039 00040 IF(n <= 0)RETURN 00041 IF (da == 0.0D0) RETURN 00042 IF(incx == 1.AND.incy == 1)GO TO 20 00043 00044 ! code for unequal increments or equal increments 00045 ! not equal to 1 00046 00047 ix = 1 00048 iy = 1 00049 IF(incx < 0)ix = (-n+1)*incx + 1 00050 IF(incy < 0)iy = (-n+1)*incy + 1 00051 DO i = 1,n 00052 dy(iy) = dy(iy) + da*dx(ix) 00053 ix = ix + incx 00054 iy = iy + incy 00055 END DO 00056 RETURN 00057 00058 ! code for both increments equal to 1 00059 00060 00061 ! clean-up loop 00062 00063 20 m = MOD(n,4) 00064 IF( m == 0 ) GO TO 40 00065 DO i = 1,m 00066 dy(i) = dy(i) + da*dx(i) 00067 END DO 00068 IF( n < 4 ) RETURN 00069 40 mp1 = m + 1 00070 DO i = mp1,n,4 00071 dy(i) = dy(i) + da*dx(i) 00072 dy(i + 1) = dy(i + 1) + da*dx(i + 1) 00073 dy(i + 2) = dy(i + 2) + da*dx(i + 2) 00074 dy(i + 3) = dy(i + 3) + da*dx(i + 3) 00075 END DO 00076 00077 END SUBROUTINE daxpy 00078 00079 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00080 00091 00092 SUBROUTINE dcopy(n,dx,incx,dy,incy) 00093 00094 DOUBLE PRECISION :: dx(*),dy(*) 00095 INTEGER :: i,incx,incy,ix,iy,m,mp1,n 00096 00097 IF(n <= 0)RETURN 00098 IF(incx == 1.AND.incy == 1)GO TO 20 00099 00100 ! code for unequal increments or equal increments 00101 ! not equal to 1 00102 00103 ix = 1 00104 iy = 1 00105 IF(incx < 0)ix = (-n+1)*incx + 1 00106 IF(incy < 0)iy = (-n+1)*incy + 1 00107 DO i = 1,n 00108 dy(iy) = dx(ix) 00109 ix = ix + incx 00110 iy = iy + incy 00111 END DO 00112 RETURN 00113 00114 ! code for both increments equal to 1 00115 00116 00117 ! clean-up loop 00118 00119 20 m = MOD(n,7) 00120 IF( m == 0 ) GO TO 40 00121 DO i = 1,m 00122 dy(i) = dx(i) 00123 END DO 00124 IF( n < 7 ) RETURN 00125 40 mp1 = m + 1 00126 DO i = mp1,n,7 00127 dy(i) = dx(i) 00128 dy(i + 1) = dx(i + 1) 00129 dy(i + 2) = dx(i + 2) 00130 dy(i + 3) = dx(i + 3) 00131 dy(i + 4) = dx(i + 4) 00132 dy(i + 5) = dx(i + 5) 00133 dy(i + 6) = dx(i + 6) 00134 END DO 00135 00136 END SUBROUTINE dcopy 00137 00138 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00139 00151 00152 DOUBLE PRECISION FUNCTION ddot(n,dx,incx,dy,incy) 00153 00154 DOUBLE PRECISION :: dx(*),dy(*),dtemp 00155 INTEGER :: i,incx,incy,ix,iy,m,mp1,n 00156 00157 ddot = 0.0D0 00158 dtemp = 0.0D0 00159 IF(n <= 0)RETURN 00160 IF(incx == 1.AND.incy == 1)GO TO 20 00161 00162 ! code for unequal increments or equal increments 00163 ! not equal to 1 00164 00165 ix = 1 00166 iy = 1 00167 IF(incx < 0)ix = (-n+1)*incx + 1 00168 IF(incy < 0)iy = (-n+1)*incy + 1 00169 DO i = 1,n 00170 dtemp = dtemp + dx(ix)*dy(iy) 00171 ix = ix + incx 00172 iy = iy + incy 00173 END DO 00174 ddot = dtemp 00175 RETURN 00176 00177 ! code for both increments equal to 1 00178 00179 00180 ! clean-up loop 00181 00182 20 m = MOD(n,5) 00183 IF( m == 0 ) GO TO 40 00184 DO i = 1,m 00185 dtemp = dtemp + dx(i)*dy(i) 00186 END DO 00187 IF( n < 5 ) GO TO 60 00188 40 mp1 = m + 1 00189 DO i = mp1,n,5 00190 dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + & 00191 dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) 00192 END DO 00193 60 ddot = dtemp 00194 00195 END FUNCTION ddot 00196 00197 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00198 00211 00212 DOUBLE PRECISION FUNCTION dnrm2 ( n, x, incx ) 00213 00214 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 00215 INTEGER :: incx, n 00216 DOUBLE PRECISION :: x(*) 00217 00218 !!! double precision s1flmx 00219 PARAMETER (one = 1.0D+0, zero = 0.0D+0 ) 00220 DOUBLE PRECISION :: norm 00221 INTRINSIC ABS 00222 ! ------------------------------------------------------------------ 00223 ! flmax = s1flmx( ) 00224 flmax = 1.0D+50 00225 00226 IF ( n < 1) THEN 00227 norm = zero 00228 00229 ELSE IF (n == 1) THEN 00230 norm = ABS( x(1) ) 00231 00232 ELSE 00233 scale = zero 00234 ssq = one 00235 00236 DO ix = 1, 1+(n-1)*incx, incx 00237 00238 IF (x(ix) /= zero) THEN 00239 absxi = ABS( x(ix) ) 00240 00241 IF (scale < absxi) THEN 00242 ssq = one + ssq*(scale/absxi)**2 00243 scale = absxi 00244 ELSE 00245 ssq = ssq + (absxi/scale)**2 00246 END IF 00247 END IF 00248 END DO 00249 00250 sqt = SQRT( ssq ) 00251 IF (scale < flmax/sqt) THEN 00252 norm = scale*sqt 00253 ELSE 00254 norm = flmax 00255 END IF 00256 END IF 00257 00258 dnrm2 = norm 00259 00260 END FUNCTION dnrm2 00261 00262 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00263 00273 00274 SUBROUTINE daxpy2( n, a, x, y, z ) 00275 00276 IMPLICIT NONE 00277 INTEGER :: n 00278 DOUBLE PRECISION :: a, x(n), y(n), z(n) 00279 INTEGER :: i 00280 00281 DO i = 1, n 00282 z(i) = a*x(i) + y(i) 00283 END DO 00284 00285 END SUBROUTINE daxpy2 00286 00287 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00288 00296 00297 SUBROUTINE dload2( n, const, x ) 00298 00299 IMPLICIT NONE 00300 INTEGER :: n 00301 DOUBLE PRECISION :: const, x(n) 00302 00303 ! ------------------------------------------------------------------ 00304 ! dload2 00305 ! ------------------------------------------------------------------ 00306 00307 INTEGER :: i 00308 00309 DO i = 1, n 00310 x(i) = const 00311 END DO 00312 00313 END SUBROUTINE dload2 00314 00315 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00316 00323 00324 SUBROUTINE dscal2( n, a, x, y ) 00325 00326 IMPLICIT NONE 00327 INTEGER :: n 00328 DOUBLE PRECISION :: a, x(n), y(n) 00329 00330 INTEGER :: i 00331 00332 DO i = 1, n 00333 y(i) = a*x(i) 00334 END DO 00335 00336 END SUBROUTINE dscal2