1;;; Compiled by f2cl version: 2;;; ("f2cl1.l,v 46c1f6a93b0d 2012/05/03 04:40:28 toy $" 3;;; "f2cl2.l,v 96616d88fb7e 2008/02/22 22:19:34 rtoy $" 4;;; "f2cl3.l,v 96616d88fb7e 2008/02/22 22:19:34 rtoy $" 5;;; "f2cl4.l,v 96616d88fb7e 2008/02/22 22:19:34 rtoy $" 6;;; "f2cl5.l,v 46c1f6a93b0d 2012/05/03 04:40:28 toy $" 7;;; "f2cl6.l,v 1d5cbacbb977 2008/08/24 00:56:27 rtoy $" 8;;; "macros.l,v fceac530ef0c 2011/11/26 04:02:26 toy $") 9 10;;; Using Lisp CMU Common Lisp snapshot-2012-04 (20C Unicode) 11;;; 12;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t) 13;;; (:coerce-assigns :as-needed) (:array-type ':array) 14;;; (:array-slicing t) (:declare-common nil) 15;;; (:float-format double-float)) 16 17(in-package :slatec) 18 19 20(let ((wg 21 (make-array 13 22 :element-type 'double-float 23 :initial-contents '(0.011393798501026288 24 0.026354986615032137 25 0.040939156701306316 26 0.054904695975835194 0.06803833381235691 27 0.08014070033500102 0.09102826198296365 28 0.10053594906705064 0.10851962447426365 29 0.11485825914571164 0.11945576353578477 30 0.12224244299031004 31 0.12317605372671545))) 32 (xgk 33 (make-array 26 34 :element-type 'double-float 35 :initial-contents '(0.9992621049926098 0.9955569697904981 36 0.9880357945340772 0.9766639214595175 37 0.9616149864258425 0.9429745712289743 38 0.9207471152817016 0.8949919978782753 39 0.8658470652932756 0.833442628760834 40 0.7978737979985001 0.7592592630373576 41 0.7177664068130843 0.6735663684734684 42 0.6268100990103174 0.577662930241223 43 0.5263252843347191 0.473002731445715 44 0.4178853821930377 0.36117230580938786 45 0.30308953893110785 0.24386688372098844 46 0.1837189394210489 0.1228646926107104 47 0.06154448300568508 0.0))) 48 (wgk 49 (make-array 26 50 :element-type 'double-float 51 :initial-contents '(0.001987383892330316 52 0.005561932135356714 53 0.009473973386174152 54 0.013236229195571676 0.0168478177091283 55 0.020435371145882834 56 0.024009945606953215 0.02747531758785174 57 0.030792300167387487 58 0.034002130274329335 0.03711627148341554 59 0.04008382550403238 0.04287284502017005 60 0.04550291304992179 0.04798253713883671 61 0.05027767908071567 0.05236288580640747 62 0.05425112988854549 0.055950811220412316 63 0.057437116361567835 64 0.058689680022394206 0.05972034032417406 65 0.06053945537604586 0.061128509717053046 66 0.061471189871425316 67 0.061580818067832936)))) 68 (declare (type (array double-float (13)) wg) 69 (type (array double-float (26)) xgk wgk)) 70 (defun dqk51 (f a b result abserr resabs resasc) 71 (declare (type (double-float) resasc resabs abserr result b a)) 72 (prog ((fv1 (make-array 25 :element-type 'double-float)) 73 (fv2 (make-array 25 :element-type 'double-float)) (j 0) (jtw 0) 74 (jtwm1 0) (absc 0.0) (centr 0.0) (dhlgth 0.0) (epmach 0.0) (fc 0.0) 75 (fsum 0.0) (fval1 0.0) (fval2 0.0) (hlgth 0.0) (resg 0.0) (resk 0.0) 76 (reskh 0.0) (uflow 0.0)) 77 (declare (type (array double-float (25)) fv2 fv1) 78 (type (double-float) uflow reskh resk resg hlgth fval2 fval1 79 fsum fc epmach dhlgth centr absc) 80 (type (f2cl-lib:integer4) jtwm1 jtw j)) 81 (setf epmach (f2cl-lib:d1mach 4)) 82 (setf uflow (f2cl-lib:d1mach 1)) 83 (setf centr (* 0.5 (+ a b))) 84 (setf hlgth (* 0.5 (- b a))) 85 (setf dhlgth (abs hlgth)) 86 (setf fc 87 (multiple-value-bind (ret-val var-0) 88 (funcall f centr) 89 (declare (ignore)) 90 (when var-0 91 (setf centr var-0)) 92 ret-val)) 93 (setf resg (* (f2cl-lib:fref wg (13) ((1 13))) fc)) 94 (setf resk (* (f2cl-lib:fref wgk (26) ((1 26))) fc)) 95 (setf resabs (abs resk)) 96 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1)) 97 ((> j 12) nil) 98 (tagbody 99 (setf jtw (f2cl-lib:int-mul j 2)) 100 (setf absc (* hlgth (f2cl-lib:fref xgk (jtw) ((1 26))))) 101 (setf fval1 (funcall f (- centr absc))) 102 (setf fval2 (funcall f (+ centr absc))) 103 (setf (f2cl-lib:fref fv1 (jtw) ((1 25))) fval1) 104 (setf (f2cl-lib:fref fv2 (jtw) ((1 25))) fval2) 105 (setf fsum (+ fval1 fval2)) 106 (setf resg (+ resg (* (f2cl-lib:fref wg (j) ((1 13))) fsum))) 107 (setf resk (+ resk (* (f2cl-lib:fref wgk (jtw) ((1 26))) fsum))) 108 (setf resabs 109 (+ resabs 110 (* (f2cl-lib:fref wgk (jtw) ((1 26))) 111 (+ (abs fval1) (abs fval2))))) 112 label10)) 113 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1)) 114 ((> j 13) nil) 115 (tagbody 116 (setf jtwm1 (f2cl-lib:int-sub (f2cl-lib:int-mul j 2) 1)) 117 (setf absc (* hlgth (f2cl-lib:fref xgk (jtwm1) ((1 26))))) 118 (setf fval1 (funcall f (- centr absc))) 119 (setf fval2 (funcall f (+ centr absc))) 120 (setf (f2cl-lib:fref fv1 (jtwm1) ((1 25))) fval1) 121 (setf (f2cl-lib:fref fv2 (jtwm1) ((1 25))) fval2) 122 (setf fsum (+ fval1 fval2)) 123 (setf resk (+ resk (* (f2cl-lib:fref wgk (jtwm1) ((1 26))) fsum))) 124 (setf resabs 125 (+ resabs 126 (* (f2cl-lib:fref wgk (jtwm1) ((1 26))) 127 (+ (abs fval1) (abs fval2))))) 128 label15)) 129 (setf reskh (* resk 0.5)) 130 (setf resasc (* (f2cl-lib:fref wgk (26) ((1 26))) (abs (- fc reskh)))) 131 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1)) 132 ((> j 25) nil) 133 (tagbody 134 (setf resasc 135 (+ resasc 136 (* (f2cl-lib:fref wgk (j) ((1 26))) 137 (+ (abs (- (f2cl-lib:fref fv1 (j) ((1 25))) reskh)) 138 (abs (- (f2cl-lib:fref fv2 (j) ((1 25))) reskh)))))) 139 label20)) 140 (setf result (* resk hlgth)) 141 (setf resabs (* resabs dhlgth)) 142 (setf resasc (* resasc dhlgth)) 143 (setf abserr (abs (* (- resk resg) hlgth))) 144 (if (and (/= resasc 0.0) (/= abserr 0.0)) 145 (setf abserr 146 (* resasc (min 1.0 (expt (/ (* 200.0 abserr) resasc) 1.5))))) 147 (if (> resabs (/ uflow (* 50.0 epmach))) 148 (setf abserr (max (* epmach 50.0 resabs) abserr))) 149 (go end_label) 150 end_label 151 (return (values nil nil nil result abserr resabs resasc))))) 152 153(in-package #-gcl #:cl-user #+gcl "CL-USER") 154#+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or)) 155(eval-when (:load-toplevel :compile-toplevel :execute) 156 (setf (gethash 'fortran-to-lisp::dqk51 fortran-to-lisp::*f2cl-function-info*) 157 (fortran-to-lisp::make-f2cl-finfo 158 :arg-types '(t (double-float) (double-float) (double-float) 159 (double-float) (double-float) (double-float)) 160 :return-values '(nil nil nil fortran-to-lisp::result 161 fortran-to-lisp::abserr fortran-to-lisp::resabs 162 fortran-to-lisp::resasc) 163 :calls '(fortran-to-lisp::d1mach)))) 164 165