share/hompack/lisp/tangns.lisp

   1 ;;; Compiled by f2cl version:
   2 ;;; ("f2cl1.l,v 95098eb54f13 2013/04/01 00:45:16 toy $"
   3 ;;;  "f2cl2.l,v 95098eb54f13 2013/04/01 00:45:16 toy $"
   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 95098eb54f13 2013/04/01 00:45:16 toy $"
   7 ;;;  "f2cl6.l,v 1d5cbacbb977 2008/08/24 00:56:27 rtoy $"
   8 ;;;  "macros.l,v 1409c1352feb 2013/03/24 20:44:50 toy $")
   9
  10 ;;; Using Lisp CMU Common Lisp snapshot-2020-04 (21D 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 "HOMPACK")
  18
  19
  20 (defun tangns
  21        (rholen y yp tz ypold a qr lenqr pivot pp rhovec work nfe n iflag par
  22         ipar)
  23   (declare (type (array f2cl-lib:integer4 (*)) ipar)
  24            (type (array double-float (*)) par)
  25            (type (array f2cl-lib:integer4 (*)) pivot)
  26            (type (f2cl-lib:integer4) iflag n nfe lenqr)
  27            (type (array double-float (*)) work rhovec pp qr a ypold tz yp y)
  28            (type (double-float) rholen))
  29   (f2cl-lib:with-multi-array-data
  30       ((y double-float y-%data% y-%offset%)
  31        (yp double-float yp-%data% yp-%offset%)
  32        (tz double-float tz-%data% tz-%offset%)
  33        (ypold double-float ypold-%data% ypold-%offset%)
  34        (a double-float a-%data% a-%offset%)
  35        (qr double-float qr-%data% qr-%offset%)
  36        (pp double-float pp-%data% pp-%offset%)
  37        (rhovec double-float rhovec-%data% rhovec-%offset%)
  38        (work double-float work-%data% work-%offset%)
  39        (pivot f2cl-lib:integer4 pivot-%data% pivot-%offset%)
  40        (par double-float par-%data% par-%offset%)
  41        (ipar f2cl-lib:integer4 ipar-%data% ipar-%offset%))
  42     (prog ((j 0) (np1 0) (np2 0) (n2p3 0) (n3p4 0) (n4p5 0) (sigma 0.0)
  43            (ypnorm 0.0) (lambda$ 0.0))
  44       (declare (type (double-float) lambda$ ypnorm sigma)
  45                (type (f2cl-lib:integer4) n4p5 n3p4 n2p3 np2 np1 j))
  46       (setf np1 (f2cl-lib:int-add n 1))
  47       (setf np2 (f2cl-lib:int-add n 2))
  48       (setf n2p3 (f2cl-lib:int-add (f2cl-lib:int-mul 2 n) 3))
  49       (setf n3p4 (f2cl-lib:int-add (f2cl-lib:int-mul 3 n) 4))
  50       (setf n4p5 (f2cl-lib:int-add (f2cl-lib:int-mul 4 n) 5))
  51       (setf nfe (f2cl-lib:int-add nfe 1))
  52       (setf lambda$
  53               (f2cl-lib:fref y-%data%
  54                              (np1)
  55                              ((1 (f2cl-lib:int-add n 1)))
  56                              y-%offset%))
  57       (cond
  58         ((= iflag (f2cl-lib:int-sub 2))
  59          (rhojs a lambda$ y qr lenqr pivot pp par ipar)
  60          (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
  61              (rho a lambda$ y rhovec par ipar)
  62            (declare (ignore var-0 var-2 var-3 var-4 var-5))
  63            (setf lambda$ var-1)))
  64         (t
  65          (f y pp)
  66          (dcopy n y 1 rhovec 1)
  67          (daxpy n -1.0 a 1 rhovec 1)
  68          (cond
  69            ((= iflag 0)
  70             (daxpy n -1.0 a 1 pp 1)
  71             (fjacs y qr lenqr pivot)
  72             (dscal lenqr (- lambda$) qr 1)
  73             (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
  74                           ((> j n) nil)
  75               (tagbody
  76                 (setf (f2cl-lib:fref qr-%data%
  77                                      ((f2cl-lib:fref pivot
  78                                                      (j)
  79                                                      ((1
  80                                                        (f2cl-lib:int-add n
  81                                                                          2)))))
  82                                      ((1 lenqr))
  83                                      qr-%offset%)
  84                         (+
  85                          (f2cl-lib:fref qr-%data%
  86                                         ((f2cl-lib:fref pivot
  87                                                         (j)
  88                                                         ((1
  89                                                           (f2cl-lib:int-add n
  90                                                                             2)))))
  91                                         ((1 lenqr))
  92                                         qr-%offset%)
  93                          1.0f0))
  94                label120))
  95             (daxpy n (- lambda$) pp 1 rhovec 1))
  96            (t
  97             (dscal n -1.0 pp 1)
  98             (daxpy n 1.0 rhovec 1 pp 1)
  99             (fjacs y qr lenqr pivot)
 100             (dscal lenqr lambda$ qr 1)
 101             (setf sigma (- 1.0f0 lambda$))
 102             (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 103                           ((> j n) nil)
 104               (tagbody
 105                 (setf (f2cl-lib:fref qr-%data%
 106                                      ((f2cl-lib:fref pivot
 107                                                      (j)
 108                                                      ((1
 109                                                        (f2cl-lib:int-add n
 110                                                                          2)))))
 111                                      ((1 lenqr))
 112                                      qr-%offset%)
 113                         (+
 114                          (f2cl-lib:fref qr-%data%
 115                                         ((f2cl-lib:fref pivot
 116                                                         (j)
 117                                                         ((1
 118                                                           (f2cl-lib:int-add n
 119                                                                             2)))))
 120                                         ((1 lenqr))
 121                                         qr-%offset%)
 122                          sigma))
 123                label170))
 124             (daxpy n (- lambda$) pp 1 rhovec 1)))))
 125       (if (< rholen 0.0f0) (setf rholen (dnrm2 n rhovec 1)))
 126       (dcopy (f2cl-lib:int-mul 2 np1) work 1
 127        (f2cl-lib:array-slice work-%data%
 128                              double-float
 129                              (n3p4)
 130                              ((1
 131                                (f2cl-lib:int-add
 132                                 (f2cl-lib:int-mul 8 (f2cl-lib:int-add n 1))
 133                                 lenqr)))
 134                              work-%offset%)
 135        1)
 136       (dcopy np1 ypold 1 yp 1)
 137       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7)
 138           (pcgds n qr lenqr pivot pp yp
 139            (f2cl-lib:array-slice work-%data%
 140                                  double-float
 141                                  (n2p3)
 142                                  ((1
 143                                    (f2cl-lib:int-add
 144                                     (f2cl-lib:int-mul 8 (f2cl-lib:int-add n 1))
 145                                     lenqr)))
 146                                  work-%offset%)
 147            iflag)
 148         (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6))
 149         (setf iflag var-7))
 150       (if (> iflag 0) (go end_label))
 151       (dcopy (f2cl-lib:int-mul 2 np1)
 152        (f2cl-lib:array-slice work-%data%
 153                              double-float
 154                              (n3p4)
 155                              ((1
 156                                (f2cl-lib:int-add
 157                                 (f2cl-lib:int-mul 8 (f2cl-lib:int-add n 1))
 158                                 lenqr)))
 159                              work-%offset%)
 160        1 work 1)
 161       (setf ypnorm (dnrm2 np1 yp 1))
 162       (dscal np1 (/ 1.0f0 ypnorm) yp 1)
 163       (if (< (ddot np1 yp 1 ypold 1) 0.0f0) (dscal np1 -1.0 yp 1))
 164       (dscal (f2cl-lib:int-mul 2 np1) 0.0
 165        (f2cl-lib:array-slice work-%data%
 166                              double-float
 167                              (n3p4)
 168                              ((1
 169                                (f2cl-lib:int-add
 170                                 (f2cl-lib:int-mul 8 (f2cl-lib:int-add n 1))
 171                                 lenqr)))
 172                              work-%offset%)
 173        1)
 174       (dcopy np1 ypold 1 tz 1)
 175       (multiple-value-bind
 176             (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8)
 177           (pcgns n qr lenqr pivot pp rhovec tz
 178            (f2cl-lib:array-slice work-%data%
 179                                  double-float
 180                                  (n2p3)
 181                                  ((1
 182                                    (f2cl-lib:int-add
 183                                     (f2cl-lib:int-mul 8 (f2cl-lib:int-add n 1))
 184                                     lenqr)))
 185                                  work-%offset%)
 186            iflag)
 187         (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7))
 188         (setf iflag var-8))
 189       (if (> iflag 0) (go end_label))
 190       (setf sigma (ddot np1 tz 1 yp 1))
 191       (daxpy np1 (- sigma) yp 1 tz 1)
 192       (go end_label)
 193      end_label
 194       (return
 195        (values rholen
 196                nil
 197                nil
 198                nil
 199                nil
 200                nil
 201                nil
 202                nil
 203                nil
 204                nil
 205                nil
 206                nil
 207                nfe
 208                nil
 209                iflag
 210                nil
 211                nil)))))
 212
 213 (in-package #-gcl #:cl-user #+gcl "CL-USER")
 214 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 215 (eval-when (:load-toplevel :compile-toplevel :execute)
 216   (setf (gethash 'fortran-to-lisp::tangns
 217                  fortran-to-lisp::*f2cl-function-info*)
 218           (fortran-to-lisp::make-f2cl-finfo
 219            :arg-types '((double-float) (array double-float (*))
 220                         (array double-float (*)) (array double-float (*))
 221                         (array double-float (*)) (array double-float (*))
 222                         (array double-float (*)) (fortran-to-lisp::integer4)
 223                         (array fortran-to-lisp::integer4 (*))
 224                         (array double-float (*)) (array double-float (*))
 225                         (array double-float (*)) (fortran-to-lisp::integer4)
 226                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 227                         (array double-float (*))
 228                         (array fortran-to-lisp::integer4 (*)))
 229            :return-values '(fortran-to-lisp::rholen nil nil nil nil nil nil nil
 230                             nil nil nil nil fortran-to-lisp::nfe nil
 231                             fortran-to-lisp::iflag nil nil)
 232            :calls '(fortran-to-lisp::fjacs fortran-to-lisp::f
 233                     fortran-to-lisp::rhojs fortran-to-lisp::ddot
 234                     fortran-to-lisp::dnrm2 fortran-to-lisp::dscal
 235                     fortran-to-lisp::daxpy fortran-to-lisp::dcopy
 236                     fortran-to-lisp::pcgns fortran-to-lisp::pcgds
 237                     fortran-to-lisp::rho))))
 238