share/minpack/lisp/qform.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))
  21   (declare (type (double-float) one zero))
  22   (defun qform (m n q ldq wa)
  23     (declare (type (array double-float (*)) wa q)
  24              (type (f2cl-lib:integer4) ldq n m))
  25     (f2cl-lib:with-multi-array-data
  26         ((q double-float q-%data% q-%offset%)
  27          (wa double-float wa-%data% wa-%offset%))
  28       (prog ((sum 0.0) (temp 0.0) (i 0) (j 0) (jm1 0) (k 0) (l 0) (minmn 0)
  29              (np1 0))
  30         (declare (type (f2cl-lib:integer4) np1 minmn l k jm1 j i)
  31                  (type (double-float) temp sum))
  32         '"     **********"
  33         '""
  34         '"     subroutine qform"
  35         '""
  36         '"     this subroutine proceeds from the computed qr factorization of"
  37         '"     an m by n matrix a to accumulate the m by m orthogonal matrix"
  38         '"     q from its factored form."
  39         '""
  40         '"     the subroutine statement is"
  41         '""
  42         '"       subroutine qform(m,n,q,ldq,wa)"
  43         '""
  44         '"     where"
  45         '""
  46         '"       m is a positive integer input variable set to the number"
  47         '"         of rows of a and the order of q."
  48         '""
  49         '"       n is a positive integer input variable set to the number"
  50         '"         of columns of a."
  51         '""
  52         '"       q is an m by m array. on input the full lower trapezoid in"
  53         '"         the first min(m,n) columns of q contains the factored form."
  54         '"         on output q has been accumulated into a square matrix."
  55         '""
  56         '"       ldq is a positive integer input variable not less than m"
  57         '"         which specifies the leading dimension of the array q."
  58         '""
  59         '"       wa is a work array of length m."
  60         '""
  61         '"     subprograms called"
  62         '""
  63         '"       fortran-supplied ... min0"
  64         '""
  65         '"     argonne national laboratory. minpack project. march 1980."
  66         '"     burton s. garbow, kenneth e. hillstrom, jorge j. more"
  67         '""
  68         '"     **********"
  69         '""
  70         '"     zero out upper triangle of q in the first min(m,n) columns."
  71         '""
  72         (setf minmn (f2cl-lib:min0 m n))
  73         (if (< minmn 2) (go label30))
  74         (f2cl-lib:fdo (j 2 (f2cl-lib:int-add j 1))
  75                       ((> j minmn) nil)
  76           (tagbody
  77             (setf jm1 (f2cl-lib:int-sub j 1))
  78             (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  79                           ((> i jm1) nil)
  80               (tagbody
  81                 (setf (f2cl-lib:fref q-%data% (i j) ((1 ldq) (1 m)) q-%offset%)
  82                         zero)
  83                label10))
  84            label20))
  85        label30
  86         '""
  87         '"     initialize remaining columns to those of the identity matrix."
  88         '""
  89         (setf np1 (f2cl-lib:int-add n 1))
  90         (if (< m np1) (go label60))
  91         (f2cl-lib:fdo (j np1 (f2cl-lib:int-add j 1))
  92                       ((> j m) nil)
  93           (tagbody
  94             (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  95                           ((> i m) nil)
  96               (tagbody
  97                 (setf (f2cl-lib:fref q-%data% (i j) ((1 ldq) (1 m)) q-%offset%)
  98                         zero)
  99                label40))
 100             (setf (f2cl-lib:fref q-%data% (j j) ((1 ldq) (1 m)) q-%offset%)
 101                     one)
 102            label50))
 103        label60
 104         '""
 105         '"     accumulate q from its factored form."
 106         '""
 107         (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 108                       ((> l minmn) nil)
 109           (tagbody
 110             (setf k (f2cl-lib:int-add (f2cl-lib:int-sub minmn l) 1))
 111             (f2cl-lib:fdo (i k (f2cl-lib:int-add i 1))
 112                           ((> i m) nil)
 113               (tagbody
 114                 (setf (f2cl-lib:fref wa-%data% (i) ((1 m)) wa-%offset%)
 115                         (f2cl-lib:fref q-%data%
 116                                        (i k)
 117                                        ((1 ldq) (1 m))
 118                                        q-%offset%))
 119                 (setf (f2cl-lib:fref q-%data% (i k) ((1 ldq) (1 m)) q-%offset%)
 120                         zero)
 121                label70))
 122             (setf (f2cl-lib:fref q-%data% (k k) ((1 ldq) (1 m)) q-%offset%)
 123                     one)
 124             (if (= (f2cl-lib:fref wa-%data% (k) ((1 m)) wa-%offset%) zero)
 125                 (go label110))
 126             (f2cl-lib:fdo (j k (f2cl-lib:int-add j 1))
 127                           ((> j m) nil)
 128               (tagbody
 129                 (setf sum zero)
 130                 (f2cl-lib:fdo (i k (f2cl-lib:int-add i 1))
 131                               ((> i m) nil)
 132                   (tagbody
 133                     (setf sum
 134                             (+ sum
 135                                (*
 136                                 (f2cl-lib:fref q-%data%
 137                                                (i j)
 138                                                ((1 ldq) (1 m))
 139                                                q-%offset%)
 140                                 (f2cl-lib:fref wa-%data%
 141                                                (i)
 142                                                ((1 m))
 143                                                wa-%offset%))))
 144                    label80))
 145                 (setf temp
 146                         (/ sum
 147                            (f2cl-lib:fref wa-%data% (k) ((1 m)) wa-%offset%)))
 148                 (f2cl-lib:fdo (i k (f2cl-lib:int-add i 1))
 149                               ((> i m) nil)
 150                   (tagbody
 151                     (setf (f2cl-lib:fref q-%data%
 152                                          (i j)
 153                                          ((1 ldq) (1 m))
 154                                          q-%offset%)
 155                             (-
 156                              (f2cl-lib:fref q-%data%
 157                                             (i j)
 158                                             ((1 ldq) (1 m))
 159                                             q-%offset%)
 160                              (* temp
 161                                 (f2cl-lib:fref wa-%data%
 162                                                (i)
 163                                                ((1 m))
 164                                                wa-%offset%))))
 165                    label90))
 166                label100))
 167            label110
 168            label120))
 169         (go end_label)
 170         '""
 171         '"     last card of subroutine qform."
 172         '""
 173        end_label
 174         (return (values nil nil nil nil nil))))))
 175
 176 (in-package #:cl-user)
 177 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 178 (eval-when (:load-toplevel :compile-toplevel :execute)
 179   (setf (gethash 'fortran-to-lisp::qform fortran-to-lisp::*f2cl-function-info*)
 180           (fortran-to-lisp::make-f2cl-finfo
 181            :arg-types '((fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 182                         (array double-float (*)) (fortran-to-lisp::integer4)
 183                         (array double-float (*)))
 184            :return-values '(nil nil nil nil nil)
 185            :calls 'nil)))
 186