share/minpack/lisp/qrfac.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) (p05 0.05) (zero 0.0))
  21   (declare (type (double-float) one p05 zero))
  22   (defun qrfac (m n a lda pivot ipvt lipvt rdiag acnorm wa)
  23     (declare (type (array f2cl-lib:integer4 (*)) ipvt)
  24              (type f2cl-lib:logical pivot)
  25              (type (array double-float (*)) wa acnorm rdiag a)
  26              (type (f2cl-lib:integer4) lipvt lda n m))
  27     (f2cl-lib:with-multi-array-data
  28         ((a double-float a-%data% a-%offset%)
  29          (rdiag double-float rdiag-%data% rdiag-%offset%)
  30          (acnorm double-float acnorm-%data% acnorm-%offset%)
  31          (wa double-float wa-%data% wa-%offset%)
  32          (ipvt f2cl-lib:integer4 ipvt-%data% ipvt-%offset%))
  33       (prog ((ajnorm 0.0) (epsmch 0.0) (sum 0.0) (temp 0.0) (i 0) (j 0) (jp1 0)
  34              (k 0) (kmax 0) (minmn 0))
  35         (declare (type (f2cl-lib:integer4) minmn kmax k jp1 j i)
  36                  (type (double-float) temp sum epsmch ajnorm))
  37         '"     **********"
  38         '""
  39         '"     subroutine qrfac"
  40         '""
  41         '"     this subroutine uses householder transformations with column"
  42         '"     pivoting (optional) to compute a qr factorization of the"
  43         '"     m by n matrix a. that is, qrfac determines an orthogonal"
  44         '"     matrix q, a permutation matrix p, and an upper trapezoidal"
  45         '"     matrix r with diagonal elements of nonincreasing magnitude,"
  46         '"     such that a*p = q*r. the householder transformation for"
  47         '"     column k, k = 1,2,...,min(m,n), is of the form"
  48         '""
  49         '"                           t"
  50         '"           i - (1/u(k))*u*u"
  51         '""
  52         '"     where u has zeros in the first k-1 positions. the form of"
  53         '"     this transformation and the method of pivoting first"
  54         '"     appeared in the corresponding linpack subroutine."
  55         '""
  56         '"     the subroutine statement is"
  57         '""
  58         '"       subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)"
  59         '""
  60         '"     where"
  61         '""
  62         '"       m is a positive integer input variable set to the number"
  63         '"         of rows of a."
  64         '""
  65         '"       n is a positive integer input variable set to the number"
  66         '"         of columns of a."
  67         '""
  68         '"       a is an m by n array. on input a contains the matrix for"
  69         '"         which the qr factorization is to be computed. on output"
  70         '"         the strict upper trapezoidal part of a contains the strict"
  71         '"         upper trapezoidal part of r, and the lower trapezoidal"
  72         '"         part of a contains a factored form of q (the non-trivial"
  73         '"         elements of the u vectors described above)."
  74         '""
  75         '"       lda is a positive integer input variable not less than m"
  76         '"         which specifies the leading dimension of the array a."
  77         '""
  78         '"       pivot is a logical input variable. if pivot is set true,"
  79         '"         then column pivoting is enforced. if pivot is set false,"
  80         '"         then no column pivoting is done."
  81         '""
  82         '"       ipvt is an integer output array of length lipvt. ipvt"
  83         '"         defines the permutation matrix p such that a*p = q*r."
  84         '"         column j of p is column ipvt(j) of the identity matrix."
  85         '"         if pivot is false, ipvt is not referenced."
  86         '""
  87         '"       lipvt is a positive integer input variable. if pivot is false,"
  88         '"         then lipvt may be as small as 1. if pivot is true, then"
  89         '"         lipvt must be at least n."
  90         '""
  91         '"       rdiag is an output array of length n which contains the"
  92         '"         diagonal elements of r."
  93         '""
  94         '"       acnorm is an output array of length n which contains the"
  95         '"         norms of the corresponding columns of the input matrix a."
  96         '"         if this information is not needed, then acnorm can coincide"
  97         '"         with rdiag."
  98         '""
  99         '"       wa is a work array of length n. if pivot is false, then wa"
 100         '"         can coincide with rdiag."
 101         '""
 102         '"     subprograms called"
 103         '""
 104         '"       minpack-supplied ... dpmpar,enorm"
 105         '""
 106         '"       fortran-supplied ... dmax1,dsqrt,min0"
 107         '""
 108         '"     argonne national laboratory. minpack project. march 1980."
 109         '"     burton s. garbow, kenneth e. hillstrom, jorge j. more"
 110         '""
 111         '"     **********"
 112         '""
 113         '"     epsmch is the machine precision."
 114         '""
 115         (setf epsmch (dpmpar 1))
 116         '""
 117         '"     compute the initial column norms and initialize several arrays."
 118         '""
 119         (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 120                       ((> j n) nil)
 121           (tagbody
 122             (setf (f2cl-lib:fref acnorm-%data% (j) ((1 n)) acnorm-%offset%)
 123                     (enorm m
 124                      (f2cl-lib:array-slice a
 125                                            double-float
 126                                            (1 j)
 127                                            ((1 lda) (1 n)))))
 128             (setf (f2cl-lib:fref rdiag-%data% (j) ((1 n)) rdiag-%offset%)
 129                     (f2cl-lib:fref acnorm-%data% (j) ((1 n)) acnorm-%offset%))
 130             (setf (f2cl-lib:fref wa-%data% (j) ((1 n)) wa-%offset%)
 131                     (f2cl-lib:fref rdiag-%data% (j) ((1 n)) rdiag-%offset%))
 132             (if pivot
 133                 (setf (f2cl-lib:fref ipvt-%data% (j) ((1 lipvt)) ipvt-%offset%)
 134                         j))
 135            label10))
 136         '""
 137         '"     reduce a to r with householder transformations."
 138         '""
 139         (setf minmn (f2cl-lib:min0 m n))
 140         (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 141                       ((> j minmn) nil)
 142           (tagbody
 143             (if (not pivot) (go label40))
 144             '""
 145             '"        bring the column of largest norm into the pivot position."
 146             '""
 147             (setf kmax j)
 148             (f2cl-lib:fdo (k j (f2cl-lib:int-add k 1))
 149                           ((> k n) nil)
 150               (tagbody
 151                 (if
 152                  (> (f2cl-lib:fref rdiag-%data% (k) ((1 n)) rdiag-%offset%)
 153                     (f2cl-lib:fref rdiag-%data% (kmax) ((1 n)) rdiag-%offset%))
 154                  (setf kmax k))
 155                label20))
 156             (if (= kmax j) (go label40))
 157             (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 158                           ((> i m) nil)
 159               (tagbody
 160                 (setf temp
 161                         (f2cl-lib:fref a-%data%
 162                                        (i j)
 163                                        ((1 lda) (1 n))
 164                                        a-%offset%))
 165                 (setf (f2cl-lib:fref a-%data% (i j) ((1 lda) (1 n)) a-%offset%)
 166                         (f2cl-lib:fref a-%data%
 167                                        (i kmax)
 168                                        ((1 lda) (1 n))
 169                                        a-%offset%))
 170                 (setf (f2cl-lib:fref a-%data%
 171                                      (i kmax)
 172                                      ((1 lda) (1 n))
 173                                      a-%offset%)
 174                         temp)
 175                label30))
 176             (setf (f2cl-lib:fref rdiag-%data% (kmax) ((1 n)) rdiag-%offset%)
 177                     (f2cl-lib:fref rdiag-%data% (j) ((1 n)) rdiag-%offset%))
 178             (setf (f2cl-lib:fref wa-%data% (kmax) ((1 n)) wa-%offset%)
 179                     (f2cl-lib:fref wa-%data% (j) ((1 n)) wa-%offset%))
 180             (setf k (f2cl-lib:fref ipvt-%data% (j) ((1 lipvt)) ipvt-%offset%))
 181             (setf (f2cl-lib:fref ipvt-%data% (j) ((1 lipvt)) ipvt-%offset%)
 182                     (f2cl-lib:fref ipvt-%data%
 183                                    (kmax)
 184                                    ((1 lipvt))
 185                                    ipvt-%offset%))
 186             (setf (f2cl-lib:fref ipvt-%data% (kmax) ((1 lipvt)) ipvt-%offset%)
 187                     k)
 188            label40
 189             '""
 190             '"        compute the householder transformation to reduce the"
 191             '"        j-th column of a to a multiple of the j-th unit vector."
 192             '""
 193             (setf ajnorm
 194                     (enorm (f2cl-lib:int-add (f2cl-lib:int-sub m j) 1)
 195                      (f2cl-lib:array-slice a
 196                                            double-float
 197                                            (j j)
 198                                            ((1 lda) (1 n)))))
 199             (if (= ajnorm zero) (go label100))
 200             (if
 201              (< (f2cl-lib:fref a-%data% (j j) ((1 lda) (1 n)) a-%offset%) zero)
 202              (setf ajnorm (- ajnorm)))
 203             (f2cl-lib:fdo (i j (f2cl-lib:int-add i 1))
 204                           ((> i m) nil)
 205               (tagbody
 206                 (setf (f2cl-lib:fref a-%data% (i j) ((1 lda) (1 n)) a-%offset%)
 207                         (/
 208                          (f2cl-lib:fref a-%data%
 209                                         (i j)
 210                                         ((1 lda) (1 n))
 211                                         a-%offset%)
 212                          ajnorm))
 213                label50))
 214             (setf (f2cl-lib:fref a-%data% (j j) ((1 lda) (1 n)) a-%offset%)
 215                     (+
 216                      (f2cl-lib:fref a-%data% (j j) ((1 lda) (1 n)) a-%offset%)
 217                      one))
 218             '""
 219             '"        apply the transformation to the remaining columns"
 220             '"        and update the norms."
 221             '""
 222             (setf jp1 (f2cl-lib:int-add j 1))
 223             (if (< n jp1) (go label100))
 224             (f2cl-lib:fdo (k jp1 (f2cl-lib:int-add k 1))
 225                           ((> k n) nil)
 226               (tagbody
 227                 (setf sum zero)
 228                 (f2cl-lib:fdo (i j (f2cl-lib:int-add i 1))
 229                               ((> i m) nil)
 230                   (tagbody
 231                     (setf sum
 232                             (+ sum
 233                                (*
 234                                 (f2cl-lib:fref a-%data%
 235                                                (i j)
 236                                                ((1 lda) (1 n))
 237                                                a-%offset%)
 238                                 (f2cl-lib:fref a-%data%
 239                                                (i k)
 240                                                ((1 lda) (1 n))
 241                                                a-%offset%))))
 242                    label60))
 243                 (setf temp
 244                         (/ sum
 245                            (f2cl-lib:fref a-%data%
 246                                           (j j)
 247                                           ((1 lda) (1 n))
 248                                           a-%offset%)))
 249                 (f2cl-lib:fdo (i j (f2cl-lib:int-add i 1))
 250                               ((> i m) nil)
 251                   (tagbody
 252                     (setf (f2cl-lib:fref a-%data%
 253                                          (i k)
 254                                          ((1 lda) (1 n))
 255                                          a-%offset%)
 256                             (-
 257                              (f2cl-lib:fref a-%data%
 258                                             (i k)
 259                                             ((1 lda) (1 n))
 260                                             a-%offset%)
 261                              (* temp
 262                                 (f2cl-lib:fref a-%data%
 263                                                (i j)
 264                                                ((1 lda) (1 n))
 265                                                a-%offset%))))
 266                    label70))
 267                 (if
 268                  (or (not pivot)
 269                      (= (f2cl-lib:fref rdiag-%data% (k) ((1 n)) rdiag-%offset%)
 270                         zero))
 271                  (go label80))
 272                 (setf temp
 273                         (/
 274                          (f2cl-lib:fref a-%data%
 275                                         (j k)
 276                                         ((1 lda) (1 n))
 277                                         a-%offset%)
 278                          (f2cl-lib:fref rdiag-%data%
 279                                         (k)
 280                                         ((1 n))
 281                                         rdiag-%offset%)))
 282                 (setf (f2cl-lib:fref rdiag-%data% (k) ((1 n)) rdiag-%offset%)
 283                         (*
 284                          (f2cl-lib:fref rdiag-%data%
 285                                         (k)
 286                                         ((1 n))
 287                                         rdiag-%offset%)
 288                          (f2cl-lib:dsqrt
 289                           (f2cl-lib:dmax1 zero (- one (expt temp 2))))))
 290                 (if
 291                  (>
 292                   (* p05
 293                      (expt
 294                       (/
 295                        (f2cl-lib:fref rdiag-%data% (k) ((1 n)) rdiag-%offset%)
 296                        (f2cl-lib:fref wa-%data% (k) ((1 n)) wa-%offset%))
 297                       2))
 298                   epsmch)
 299                  (go label80))
 300                 (setf (f2cl-lib:fref rdiag-%data% (k) ((1 n)) rdiag-%offset%)
 301                         (enorm (f2cl-lib:int-sub m j)
 302                          (f2cl-lib:array-slice a
 303                                                double-float
 304                                                (jp1 k)
 305                                                ((1 lda) (1 n)))))
 306                 (setf (f2cl-lib:fref wa-%data% (k) ((1 n)) wa-%offset%)
 307                         (f2cl-lib:fref rdiag-%data%
 308                                        (k)
 309                                        ((1 n))
 310                                        rdiag-%offset%))
 311                label80
 312                label90))
 313            label100
 314             (setf (f2cl-lib:fref rdiag-%data% (j) ((1 n)) rdiag-%offset%)
 315                     (- ajnorm))
 316            label110))
 317         (go end_label)
 318         '""
 319         '"     last card of subroutine qrfac."
 320         '""
 321        end_label
 322         (return (values nil nil nil nil nil nil nil nil nil nil))))))
 323
 324 (in-package #:cl-user)
 325 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 326 (eval-when (:load-toplevel :compile-toplevel :execute)
 327   (setf (gethash 'fortran-to-lisp::qrfac fortran-to-lisp::*f2cl-function-info*)
 328           (fortran-to-lisp::make-f2cl-finfo
 329            :arg-types '((fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 330                         (array double-float (*)) (fortran-to-lisp::integer4)
 331                         fortran-to-lisp::logical
 332                         (array fortran-to-lisp::integer4 (*))
 333                         (fortran-to-lisp::integer4) (array double-float (*))
 334                         (array double-float (*)) (array double-float (*)))
 335            :return-values '(nil nil nil nil nil nil nil nil nil nil)
 336            :calls '(fortran-to-lisp::enorm fortran-to-lisp::dpmpar))))
 337