share/minpack/lisp/enorm.lisp

   1 ;;; Compiled by f2cl version:
   2 ;;; ("f2cl1.l,v 1.215 2009/04/07 22:05:21 rtoy Exp $"
   3 ;;;  "f2cl2.l,v 1.37 2008/02/22 22:19:33 rtoy Exp $"
   4 ;;;  "f2cl3.l,v 1.6 2008/02/22 22:19:33 rtoy Exp $"
   5 ;;;  "f2cl4.l,v 1.7 2008/02/22 22:19:34 rtoy Exp $"
   6 ;;;  "f2cl5.l,v 1.200 2009/01/19 02:38:17 rtoy Exp $"
   7 ;;;  "f2cl6.l,v 1.48 2008/08/24 00:56:27 rtoy Exp $"
   8 ;;;  "macros.l,v 1.112 2009/01/08 12:57:19 rtoy Exp $")
   9
  10 ;;; Using Lisp CMU Common Lisp 19f (19F)
  11 ;;;
  12 ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls nil)
  13 ;;;           (:coerce-assigns :as-needed) (:array-type ':array)
  14 ;;;           (:array-slicing t) (:declare-common nil)
  15 ;;;           (:float-format double-float))
  16
  17 (in-package :minpack)
  18
  19
  20 (let ((one 1.0) (zero 0.0) (rdwarf 3.834e-20) (rgiant 1.304e19))
  21   (declare (type (double-float) one zero rdwarf rgiant))
  22   (defun enorm (n x)
  23     (declare (type (array double-float (*)) x) (type (f2cl-lib:integer4) n))
  24     (f2cl-lib:with-multi-array-data
  25         ((x double-float x-%data% x-%offset%))
  26       (prog ((agiant 0.0) (floatn 0.0) (s1 0.0) (s2 0.0) (s3 0.0) (xabs 0.0)
  27              (x1max 0.0) (x3max 0.0) (i 0) (enorm 0.0))
  28         (declare (type (f2cl-lib:integer4) i)
  29                  (type (double-float) enorm x3max x1max xabs s3 s2 s1 floatn
  30                                       agiant))
  31         '"     **********"
  32         '""
  33         '"     function enorm"
  34         '""
  35         '"     given an n-vector x, this function calculates the"
  36         '"     euclidean norm of x."
  37         '""
  38         '"     the euclidean norm is computed by accumulating the sum of"
  39         '"     squares in three different sums. the sums of squares for the"
  40         '"     small and large components are scaled so that no overflows"
  41         '"     occur. non-destructive underflows are permitted. underflows"
  42         '"     and overflows do not occur in the computation of the unscaled"
  43         '"     sum of squares for the intermediate components."
  44         '"     the definitions of small, intermediate and large components"
  45         '"     depend on two constants, rdwarf and rgiant. the main"
  46         '"     restrictions on these constants are that rdwarf**2 not"
  47         '"     underflow and rgiant**2 not overflow. the constants"
  48         '"     given here are suitable for every known computer."
  49         '""
  50         '"     the function statement is"
  51         '""
  52         '"       double precision function enorm(n,x)"
  53         '""
  54         '"     where"
  55         '""
  56         '"       n is a positive integer input variable."
  57         '""
  58         '"       x is an input array of length n."
  59         '""
  60         '"     subprograms called"
  61         '""
  62         '"       fortran-supplied ... dabs,dsqrt"
  63         '""
  64         '"     argonne national laboratory. minpack project. march 1980."
  65         '"     burton s. garbow, kenneth e. hillstrom, jorge j. more"
  66         '""
  67         '"     **********"
  68         (setf s1 zero)
  69         (setf s2 zero)
  70         (setf s3 zero)
  71         (setf x1max zero)
  72         (setf x3max zero)
  73         (setf floatn (coerce (the f2cl-lib:integer4 n) 'double-float))
  74         (setf agiant (/ rgiant floatn))
  75         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  76                       ((> i n) nil)
  77           (tagbody
  78             (setf xabs
  79                     (f2cl-lib:dabs
  80                      (f2cl-lib:fref x-%data% (i) ((1 n)) x-%offset%)))
  81             (if (and (> xabs rdwarf) (< xabs agiant)) (go label70))
  82             (if (<= xabs rdwarf) (go label30))
  83             '""
  84             '"              sum for large components."
  85             '""
  86             (if (<= xabs x1max) (go label10))
  87             (setf s1 (+ one (* s1 (expt (/ x1max xabs) 2))))
  88             (setf x1max xabs)
  89             (go label20)
  90            label10
  91             (setf s1 (+ s1 (expt (/ xabs x1max) 2)))
  92            label20
  93             (go label60)
  94            label30
  95             '""
  96             '"              sum for small components."
  97             '""
  98             (if (<= xabs x3max) (go label40))
  99             (setf s3 (+ one (* s3 (expt (/ x3max xabs) 2))))
 100             (setf x3max xabs)
 101             (go label50)
 102            label40
 103             (if (/= xabs zero) (setf s3 (+ s3 (expt (/ xabs x3max) 2))))
 104            label50
 105            label60
 106             (go label80)
 107            label70
 108             '""
 109             '"           sum for intermediate components."
 110             '""
 111             (setf s2 (+ s2 (expt xabs 2)))
 112            label80
 113            label90))
 114         '""
 115         '"     calculation of norm."
 116         '""
 117         (if (= s1 zero) (go label100))
 118         (setf enorm (* x1max (f2cl-lib:dsqrt (+ s1 (/ (/ s2 x1max) x1max)))))
 119         (go label130)
 120        label100
 121         (if (= s2 zero) (go label110))
 122         (if (>= s2 x3max)
 123             (setf enorm
 124                     (f2cl-lib:dsqrt
 125                      (* s2 (+ one (* (/ x3max s2) (* x3max s3)))))))
 126         (if (< s2 x3max)
 127             (setf enorm
 128                     (f2cl-lib:dsqrt (* x3max (+ (/ s2 x3max) (* x3max s3))))))
 129         (go label120)
 130        label110
 131         (setf enorm (* x3max (f2cl-lib:dsqrt s3)))
 132        label120
 133        label130
 134         (go end_label)
 135         '""
 136         '"     last card of function enorm."
 137         '""
 138        end_label
 139         (return (values enorm nil nil))))))
 140
 141 (in-package #:cl-user)
 142 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 143 (eval-when (:load-toplevel :compile-toplevel :execute)
 144   (setf (gethash 'fortran-to-lisp::enorm fortran-to-lisp::*f2cl-function-info*)
 145           (fortran-to-lisp::make-f2cl-finfo
 146            :arg-types '((fortran-to-lisp::integer4) (array double-float (*)))
 147            :return-values '(nil nil)
 148            :calls 'nil)))
 149