share/minpack/lisp/r1mpyq.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))
  21   (declare (type (double-float) one))
  22   (defun r1mpyq (m n a lda v w)
  23     (declare (type (array double-float (*)) w v a)
  24              (type (f2cl-lib:integer4) lda n m))
  25     (f2cl-lib:with-multi-array-data
  26         ((a double-float a-%data% a-%offset%)
  27          (v double-float v-%data% v-%offset%)
  28          (w double-float w-%data% w-%offset%))
  29       (prog ((cos 0.0) (sin 0.0) (temp 0.0) (i 0) (j 0) (nmj 0) (nm1 0))
  30         (declare (type (f2cl-lib:integer4) nm1 nmj j i)
  31                  (type (double-float) temp sin cos))
  32         '"     **********"
  33         '""
  34         '"     subroutine r1mpyq"
  35         '""
  36         '"     given an m by n matrix a, this subroutine computes a*q where"
  37         '"     q is the product of 2*(n - 1) transformations"
  38         '""
  39         '"           gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)"
  40         '""
  41         '"     and gv(i), gw(i) are givens rotations in the (i,n) plane which"
  42         '"     eliminate elements in the i-th and n-th planes, respectively."
  43         '"     q itself is not given, rather the information to recover the"
  44         '"     gv, gw rotations is supplied."
  45         '""
  46         '"     the subroutine statement is"
  47         '""
  48         '"       subroutine r1mpyq(m,n,a,lda,v,w)"
  49         '""
  50         '"     where"
  51         '""
  52         '"       m is a positive integer input variable set to the number"
  53         '"         of rows of a."
  54         '""
  55         '"       n is a positive integer input variable set to the number"
  56         '"         of columns of a."
  57         '""
  58         '"       a is an m by n array. on input a must contain the matrix"
  59         '"         to be postmultiplied by the orthogonal matrix q"
  60         '"         described above. on output a*q has replaced a."
  61         '""
  62         '"       lda is a positive integer input variable not less than m"
  63         '"         which specifies the leading dimension of the array a."
  64         '""
  65         '"       v is an input array of length n. v(i) must contain the"
  66         '"         information necessary to recover the givens rotation gv(i)"
  67         '"         described above."
  68         '""
  69         '"       w is an input array of length n. w(i) must contain the"
  70         '"         information necessary to recover the givens rotation gw(i)"
  71         '"         described above."
  72         '""
  73         '"     subroutines called"
  74         '""
  75         '"       fortran-supplied ... dabs,dsqrt"
  76         '""
  77         '"     argonne national laboratory. minpack project. march 1980."
  78         '"     burton s. garbow, kenneth e. hillstrom, jorge j. more"
  79         '""
  80         '"     **********"
  81         '""
  82         '"     apply the first set of givens rotations to a."
  83         '""
  84         (setf nm1 (f2cl-lib:int-sub n 1))
  85         (if (< nm1 1) (go label50))
  86         (f2cl-lib:fdo (nmj 1 (f2cl-lib:int-add nmj 1))
  87                       ((> nmj nm1) nil)
  88           (tagbody
  89             (setf j (f2cl-lib:int-sub n nmj))
  90             (if
  91              (> (f2cl-lib:dabs (f2cl-lib:fref v-%data% (j) ((1 n)) v-%offset%))
  92                 one)
  93              (setf cos
  94                      (/ one (f2cl-lib:fref v-%data% (j) ((1 n)) v-%offset%))))
  95             (if
  96              (> (f2cl-lib:dabs (f2cl-lib:fref v-%data% (j) ((1 n)) v-%offset%))
  97                 one)
  98              (setf sin (f2cl-lib:dsqrt (- one (expt cos 2)))))
  99             (if
 100              (<=
 101               (f2cl-lib:dabs (f2cl-lib:fref v-%data% (j) ((1 n)) v-%offset%))
 102               one)
 103              (setf sin (f2cl-lib:fref v-%data% (j) ((1 n)) v-%offset%)))
 104             (if
 105              (<=
 106               (f2cl-lib:dabs (f2cl-lib:fref v-%data% (j) ((1 n)) v-%offset%))
 107               one)
 108              (setf cos (f2cl-lib:dsqrt (- one (expt sin 2)))))
 109             (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 110                           ((> i m) nil)
 111               (tagbody
 112                 (setf temp
 113                         (-
 114                          (* cos
 115                             (f2cl-lib:fref a-%data%
 116                                            (i j)
 117                                            ((1 lda) (1 n))
 118                                            a-%offset%))
 119                          (* sin
 120                             (f2cl-lib:fref a-%data%
 121                                            (i n)
 122                                            ((1 lda) (1 n))
 123                                            a-%offset%))))
 124                 (setf (f2cl-lib:fref a-%data% (i n) ((1 lda) (1 n)) a-%offset%)
 125                         (+
 126                          (* sin
 127                             (f2cl-lib:fref a-%data%
 128                                            (i j)
 129                                            ((1 lda) (1 n))
 130                                            a-%offset%))
 131                          (* cos
 132                             (f2cl-lib:fref a-%data%
 133                                            (i n)
 134                                            ((1 lda) (1 n))
 135                                            a-%offset%))))
 136                 (setf (f2cl-lib:fref a-%data% (i j) ((1 lda) (1 n)) a-%offset%)
 137                         temp)
 138                label10))
 139            label20))
 140         '""
 141         '"     apply the second set of givens rotations to a."
 142         '""
 143         (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 144                       ((> j nm1) nil)
 145           (tagbody
 146             (if
 147              (> (f2cl-lib:dabs (f2cl-lib:fref w-%data% (j) ((1 n)) w-%offset%))
 148                 one)
 149              (setf cos
 150                      (/ one (f2cl-lib:fref w-%data% (j) ((1 n)) w-%offset%))))
 151             (if
 152              (> (f2cl-lib:dabs (f2cl-lib:fref w-%data% (j) ((1 n)) w-%offset%))
 153                 one)
 154              (setf sin (f2cl-lib:dsqrt (- one (expt cos 2)))))
 155             (if
 156              (<=
 157               (f2cl-lib:dabs (f2cl-lib:fref w-%data% (j) ((1 n)) w-%offset%))
 158               one)
 159              (setf sin (f2cl-lib:fref w-%data% (j) ((1 n)) w-%offset%)))
 160             (if
 161              (<=
 162               (f2cl-lib:dabs (f2cl-lib:fref w-%data% (j) ((1 n)) w-%offset%))
 163               one)
 164              (setf cos (f2cl-lib:dsqrt (- one (expt sin 2)))))
 165             (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 166                           ((> i m) nil)
 167               (tagbody
 168                 (setf temp
 169                         (+
 170                          (* cos
 171                             (f2cl-lib:fref a-%data%
 172                                            (i j)
 173                                            ((1 lda) (1 n))
 174                                            a-%offset%))
 175                          (* sin
 176                             (f2cl-lib:fref a-%data%
 177                                            (i n)
 178                                            ((1 lda) (1 n))
 179                                            a-%offset%))))
 180                 (setf (f2cl-lib:fref a-%data% (i n) ((1 lda) (1 n)) a-%offset%)
 181                         (+
 182                          (* (- sin)
 183                             (f2cl-lib:fref a-%data%
 184                                            (i j)
 185                                            ((1 lda) (1 n))
 186                                            a-%offset%))
 187                          (* cos
 188                             (f2cl-lib:fref a-%data%
 189                                            (i n)
 190                                            ((1 lda) (1 n))
 191                                            a-%offset%))))
 192                 (setf (f2cl-lib:fref a-%data% (i j) ((1 lda) (1 n)) a-%offset%)
 193                         temp)
 194                label30))
 195            label40))
 196        label50
 197         (go end_label)
 198         '""
 199         '"     last card of subroutine r1mpyq."
 200         '""
 201        end_label
 202         (return (values nil nil nil nil nil nil))))))
 203
 204 (in-package #:cl-user)
 205 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 206 (eval-when (:load-toplevel :compile-toplevel :execute)
 207   (setf (gethash 'fortran-to-lisp::r1mpyq
 208                  fortran-to-lisp::*f2cl-function-info*)
 209           (fortran-to-lisp::make-f2cl-finfo
 210            :arg-types '((fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 211                         (array double-float (*)) (fortran-to-lisp::integer4)
 212                         (array double-float (*)) (array double-float (*)))
 213            :return-values '(nil nil nil nil nil nil)
 214            :calls 'nil)))
 215