share/minpack/lisp/r1updt.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) (p5 0.5) (p25 0.25) (zero 0.0))
  21   (declare (type (double-float) one p5 p25 zero))
  22   (defun r1updt (m n s ls u v w sing)
  23     (declare (type f2cl-lib:logical sing)
  24              (type (array double-float (*)) w v u s)
  25              (type (f2cl-lib:integer4) ls n m))
  26     (f2cl-lib:with-multi-array-data
  27         ((s double-float s-%data% s-%offset%)
  28          (u double-float u-%data% u-%offset%)
  29          (v double-float v-%data% v-%offset%)
  30          (w double-float w-%data% w-%offset%))
  31       (prog ((cos 0.0) (cotan 0.0) (giant 0.0) (sin 0.0) (tan 0.0) (tau 0.0)
  32              (temp 0.0) (i 0) (j 0) (jj 0) (l 0) (nmj 0) (nm1 0))
  33         (declare (type (f2cl-lib:integer4) nm1 nmj l jj j i)
  34                  (type (double-float) temp tau tan sin giant cotan cos))
  35         '"     **********"
  36         '""
  37         '"     subroutine r1updt"
  38         '""
  39         '"     given an m by n lower trapezoidal matrix s, an m-vector u,"
  40         '"     and an n-vector v, the problem is to determine an"
  41         '"     orthogonal matrix q such that"
  42         '""
  43         '"                   t"
  44         '"           (s + u*v )*q"
  45         '""
  46         '"     is again lower trapezoidal."
  47         '""
  48         '"     this subroutine determines q as the product of 2*(n - 1)"
  49         '"     transformations"
  50         '""
  51         '"           gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)"
  52         '""
  53         '"     where gv(i), gw(i) are givens rotations in the (i,n) plane"
  54         '"     which eliminate elements in the i-th and n-th planes,"
  55         '"     respectively. q itself is not accumulated, rather the"
  56         '"     information to recover the gv, gw rotations is returned."
  57         '""
  58         '"     the subroutine statement is"
  59         '""
  60         '"       subroutine r1updt(m,n,s,ls,u,v,w,sing)"
  61         '""
  62         '"     where"
  63         '""
  64         '"       m is a positive integer input variable set to the number"
  65         '"         of rows of s."
  66         '""
  67         '"       n is a positive integer input variable set to the number"
  68         '"         of columns of s. n must not exceed m."
  69         '""
  70         '"       s is an array of length ls. on input s must contain the lower"
  71         '"         trapezoidal matrix s stored by columns. on output s contains"
  72         '"         the lower trapezoidal matrix produced as described above."
  73         '""
  74         '"       ls is a positive integer input variable not less than"
  75         '"         (n*(2*m-n+1))/2."
  76         '""
  77         '"       u is an input array of length m which must contain the"
  78         '"         vector u."
  79         '""
  80         '"       v is an array of length n. on input v must contain the vector"
  81         '"         v. on output v(i) contains the information necessary to"
  82         '"         recover the givens rotation gv(i) described above."
  83         '""
  84         '"       w is an output array of length m. w(i) contains information"
  85         '"         necessary to recover the givens rotation gw(i) described"
  86         '"         above."
  87         '""
  88         '"       sing is a logical output variable. sing is set true if any"
  89         '"         of the diagonal elements of the output s are zero. otherwise"
  90         '"         sing is set false."
  91         '""
  92         '"     subprograms called"
  93         '""
  94         '"       minpack-supplied ... dpmpar"
  95         '""
  96         '"       fortran-supplied ... dabs,dsqrt"
  97         '""
  98         '"     argonne national laboratory. minpack project. march 1980."
  99         '"     burton s. garbow, kenneth e. hillstrom, jorge j. more,"
 100         '"     john l. nazareth"
 101         '""
 102         '"     **********"
 103         '""
 104         '"     giant is the largest magnitude."
 105         '""
 106         (setf giant (dpmpar 3))
 107         '""
 108         '"     initialize the diagonal element pointer."
 109         '""
 110         (setf jj
 111                 (-
 112                  (the f2cl-lib:integer4 (truncate (* n (+ (- (* 2 m) n) 1)) 2))
 113                  (f2cl-lib:int-sub m n)))
 114         '""
 115         '"     move the nontrivial part of the last column of s into w."
 116         '""
 117         (setf l jj)
 118         (f2cl-lib:fdo (i n (f2cl-lib:int-add i 1))
 119                       ((> i m) nil)
 120           (tagbody
 121             (setf (f2cl-lib:fref w-%data% (i) ((1 m)) w-%offset%)
 122                     (f2cl-lib:fref s-%data% (l) ((1 ls)) s-%offset%))
 123             (setf l (f2cl-lib:int-add l 1))
 124            label10))
 125         '""
 126         '"     rotate the vector v into a multiple of the n-th unit vector"
 127         '"     in such a way that a spike is introduced into w."
 128         '""
 129         (setf nm1 (f2cl-lib:int-sub n 1))
 130         (if (< nm1 1) (go label70))
 131         (f2cl-lib:fdo (nmj 1 (f2cl-lib:int-add nmj 1))
 132                       ((> nmj nm1) nil)
 133           (tagbody
 134             (setf j (f2cl-lib:int-sub n nmj))
 135             (setf jj
 136                     (f2cl-lib:int-sub jj
 137                                       (f2cl-lib:int-add (f2cl-lib:int-sub m j)
 138                                                         1)))
 139             (setf (f2cl-lib:fref w-%data% (j) ((1 m)) w-%offset%) zero)
 140             (if (= (f2cl-lib:fref v-%data% (j) ((1 n)) v-%offset%) zero)
 141                 (go label50))
 142             '""
 143             '"        determine a givens rotation which eliminates the"
 144             '"        j-th element of v."
 145             '""
 146             (if
 147              (>=
 148               (f2cl-lib:dabs (f2cl-lib:fref v-%data% (n) ((1 n)) v-%offset%))
 149               (f2cl-lib:dabs (f2cl-lib:fref v-%data% (j) ((1 n)) v-%offset%)))
 150              (go label20))
 151             (setf cotan
 152                     (/ (f2cl-lib:fref v-%data% (n) ((1 n)) v-%offset%)
 153                        (f2cl-lib:fref v-%data% (j) ((1 n)) v-%offset%)))
 154             (setf sin (/ p5 (f2cl-lib:dsqrt (+ p25 (* p25 (expt cotan 2))))))
 155             (setf cos (* sin cotan))
 156             (setf tau one)
 157             (if (> (* (f2cl-lib:dabs cos) giant) one) (setf tau (/ one cos)))
 158             (go label30)
 159            label20
 160             (setf tan
 161                     (/ (f2cl-lib:fref v-%data% (j) ((1 n)) v-%offset%)
 162                        (f2cl-lib:fref v-%data% (n) ((1 n)) v-%offset%)))
 163             (setf cos (/ p5 (f2cl-lib:dsqrt (+ p25 (* p25 (expt tan 2))))))
 164             (setf sin (* cos tan))
 165             (setf tau sin)
 166            label30
 167             '""
 168             '"        apply the transformation to v and store the information"
 169             '"        necessary to recover the givens rotation."
 170             '""
 171             (setf (f2cl-lib:fref v-%data% (n) ((1 n)) v-%offset%)
 172                     (+ (* sin (f2cl-lib:fref v-%data% (j) ((1 n)) v-%offset%))
 173                        (* cos
 174                           (f2cl-lib:fref v-%data% (n) ((1 n)) v-%offset%))))
 175             (setf (f2cl-lib:fref v-%data% (j) ((1 n)) v-%offset%) tau)
 176             '""
 177             '"        apply the transformation to s and extend the spike in w."
 178             '""
 179             (setf l jj)
 180             (f2cl-lib:fdo (i j (f2cl-lib:int-add i 1))
 181                           ((> i m) nil)
 182               (tagbody
 183                 (setf temp
 184                         (-
 185                          (* cos
 186                             (f2cl-lib:fref s-%data% (l) ((1 ls)) s-%offset%))
 187                          (* sin
 188                             (f2cl-lib:fref w-%data% (i) ((1 m)) w-%offset%))))
 189                 (setf (f2cl-lib:fref w-%data% (i) ((1 m)) w-%offset%)
 190                         (+
 191                          (* sin
 192                             (f2cl-lib:fref s-%data% (l) ((1 ls)) s-%offset%))
 193                          (* cos
 194                             (f2cl-lib:fref w-%data% (i) ((1 m)) w-%offset%))))
 195                 (setf (f2cl-lib:fref s-%data% (l) ((1 ls)) s-%offset%) temp)
 196                 (setf l (f2cl-lib:int-add l 1))
 197                label40))
 198            label50
 199            label60))
 200        label70
 201         '""
 202         '"     add the spike from the rank 1 update to w."
 203         '""
 204         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 205                       ((> i m) nil)
 206           (tagbody
 207             (setf (f2cl-lib:fref w-%data% (i) ((1 m)) w-%offset%)
 208                     (+ (f2cl-lib:fref w-%data% (i) ((1 m)) w-%offset%)
 209                        (* (f2cl-lib:fref v-%data% (n) ((1 n)) v-%offset%)
 210                           (f2cl-lib:fref u-%data% (i) ((1 m)) u-%offset%))))
 211            label80))
 212         '""
 213         '"     eliminate the spike."
 214         '""
 215         (setf sing f2cl-lib:%false%)
 216         (if (< nm1 1) (go label140))
 217         (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 218                       ((> j nm1) nil)
 219           (tagbody
 220             (if (= (f2cl-lib:fref w-%data% (j) ((1 m)) w-%offset%) zero)
 221                 (go label120))
 222             '""
 223             '"        determine a givens rotation which eliminates the"
 224             '"        j-th element of the spike."
 225             '""
 226             (if
 227              (>=
 228               (f2cl-lib:dabs (f2cl-lib:fref s-%data% (jj) ((1 ls)) s-%offset%))
 229               (f2cl-lib:dabs (f2cl-lib:fref w-%data% (j) ((1 m)) w-%offset%)))
 230              (go label90))
 231             (setf cotan
 232                     (/ (f2cl-lib:fref s-%data% (jj) ((1 ls)) s-%offset%)
 233                        (f2cl-lib:fref w-%data% (j) ((1 m)) w-%offset%)))
 234             (setf sin (/ p5 (f2cl-lib:dsqrt (+ p25 (* p25 (expt cotan 2))))))
 235             (setf cos (* sin cotan))
 236             (setf tau one)
 237             (if (> (* (f2cl-lib:dabs cos) giant) one) (setf tau (/ one cos)))
 238             (go label100)
 239            label90
 240             (setf tan
 241                     (/ (f2cl-lib:fref w-%data% (j) ((1 m)) w-%offset%)
 242                        (f2cl-lib:fref s-%data% (jj) ((1 ls)) s-%offset%)))
 243             (setf cos (/ p5 (f2cl-lib:dsqrt (+ p25 (* p25 (expt tan 2))))))
 244             (setf sin (* cos tan))
 245             (setf tau sin)
 246            label100
 247             '""
 248             '"        apply the transformation to s and reduce the spike in w."
 249             '""
 250             (setf l jj)
 251             (f2cl-lib:fdo (i j (f2cl-lib:int-add i 1))
 252                           ((> i m) nil)
 253               (tagbody
 254                 (setf temp
 255                         (+
 256                          (* cos
 257                             (f2cl-lib:fref s-%data% (l) ((1 ls)) s-%offset%))
 258                          (* sin
 259                             (f2cl-lib:fref w-%data% (i) ((1 m)) w-%offset%))))
 260                 (setf (f2cl-lib:fref w-%data% (i) ((1 m)) w-%offset%)
 261                         (+
 262                          (* (- sin)
 263                             (f2cl-lib:fref s-%data% (l) ((1 ls)) s-%offset%))
 264                          (* cos
 265                             (f2cl-lib:fref w-%data% (i) ((1 m)) w-%offset%))))
 266                 (setf (f2cl-lib:fref s-%data% (l) ((1 ls)) s-%offset%) temp)
 267                 (setf l (f2cl-lib:int-add l 1))
 268                label110))
 269             '""
 270             '"        store the information necessary to recover the"
 271             '"        givens rotation."
 272             '""
 273             (setf (f2cl-lib:fref w-%data% (j) ((1 m)) w-%offset%) tau)
 274            label120
 275             '""
 276             '"        test for zero diagonal elements in the output s."
 277             '""
 278             (if (= (f2cl-lib:fref s-%data% (jj) ((1 ls)) s-%offset%) zero)
 279                 (setf sing f2cl-lib:%true%))
 280             (setf jj
 281                     (f2cl-lib:int-add jj
 282                                       (f2cl-lib:int-add (f2cl-lib:int-sub m j)
 283                                                         1)))
 284            label130))
 285        label140
 286         '""
 287         '"     move w back into the last column of the output s."
 288         '""
 289         (setf l jj)
 290         (f2cl-lib:fdo (i n (f2cl-lib:int-add i 1))
 291                       ((> i m) nil)
 292           (tagbody
 293             (setf (f2cl-lib:fref s-%data% (l) ((1 ls)) s-%offset%)
 294                     (f2cl-lib:fref w-%data% (i) ((1 m)) w-%offset%))
 295             (setf l (f2cl-lib:int-add l 1))
 296            label150))
 297         (if (= (f2cl-lib:fref s-%data% (jj) ((1 ls)) s-%offset%) zero)
 298             (setf sing f2cl-lib:%true%))
 299         (go end_label)
 300         '""
 301         '"     last card of subroutine r1updt."
 302         '""
 303        end_label
 304         (return (values nil nil nil nil nil nil nil sing))))))
 305
 306 (in-package #:cl-user)
 307 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 308 (eval-when (:load-toplevel :compile-toplevel :execute)
 309   (setf (gethash 'fortran-to-lisp::r1updt
 310                  fortran-to-lisp::*f2cl-function-info*)
 311           (fortran-to-lisp::make-f2cl-finfo
 312            :arg-types '((fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 313                         (array double-float (*)) (fortran-to-lisp::integer4)
 314                         (array double-float (*)) (array double-float (*))
 315                         (array double-float (*)) fortran-to-lisp::logical)
 316            :return-values '(nil nil nil nil nil nil nil fortran-to-lisp::sing)
 317            :calls '(fortran-to-lisp::dpmpar))))
 318