share/hompack/lisp/dcpose.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 dcpose (ndim n qr alpha pivot ierr y sum)
  21   (declare (type (array double-float (*)) sum y)
  22            (type (array f2cl-lib:integer4 (*)) pivot)
  23            (type (array double-float (*)) alpha qr)
  24            (type (f2cl-lib:integer4) ierr n ndim))
  25   (f2cl-lib:with-multi-array-data
  26       ((qr double-float qr-%data% qr-%offset%)
  27        (alpha double-float alpha-%data% alpha-%offset%)
  28        (pivot f2cl-lib:integer4 pivot-%data% pivot-%offset%)
  29        (y double-float y-%data% y-%offset%)
  30        (sum double-float sum-%data% sum-%offset%))
  31     (prog ((beta 0.0) (sigma 0.0) (alphak 0.0) (qrkk 0.0) (i 0) (j 0) (jbar 0)
  32            (k 0) (kp1 0) (np1 0))
  33       (declare (type (f2cl-lib:integer4) np1 kp1 k jbar j i)
  34                (type (double-float) qrkk alphak sigma beta))
  35       (setf ierr 0)
  36       (setf np1 (f2cl-lib:int-add n 1))
  37       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
  38                     ((> j np1) nil)
  39         (tagbody
  40           (setf (f2cl-lib:fref sum-%data% (j) ((1 1)) sum-%offset%)
  41                   (ddot n
  42                    (f2cl-lib:array-slice qr-%data%
  43                                          double-float
  44                                          (1 j)
  45                                          ((1 ndim) (1 1))
  46                                          qr-%offset%)
  47                    1
  48                    (f2cl-lib:array-slice qr-%data%
  49                                          double-float
  50                                          (1 j)
  51                                          ((1 ndim) (1 1))
  52                                          qr-%offset%)
  53                    1))
  54          label20
  55           (setf (f2cl-lib:fref pivot-%data% (j) ((1 1)) pivot-%offset%) j)))
  56       (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
  57                     ((> k n) nil)
  58         (tagbody
  59           (setf sigma (f2cl-lib:fref sum-%data% (k) ((1 1)) sum-%offset%))
  60           (setf jbar k)
  61           (setf kp1 (f2cl-lib:int-add k 1))
  62           (f2cl-lib:fdo (j kp1 (f2cl-lib:int-add j 1))
  63                         ((> j np1) nil)
  64             (tagbody
  65               (if
  66                (>= sigma (f2cl-lib:fref sum-%data% (j) ((1 1)) sum-%offset%))
  67                (go label40))
  68               (setf sigma (f2cl-lib:fref sum-%data% (j) ((1 1)) sum-%offset%))
  69               (setf jbar j)
  70              label40))
  71           (if (= jbar k) (go label70))
  72           (setf i (f2cl-lib:fref pivot-%data% (k) ((1 1)) pivot-%offset%))
  73           (setf (f2cl-lib:fref pivot-%data% (k) ((1 1)) pivot-%offset%)
  74                   (f2cl-lib:fref pivot-%data% (jbar) ((1 1)) pivot-%offset%))
  75           (setf (f2cl-lib:fref pivot-%data% (jbar) ((1 1)) pivot-%offset%) i)
  76           (setf (f2cl-lib:fref sum-%data% (jbar) ((1 1)) sum-%offset%)
  77                   (f2cl-lib:fref sum-%data% (k) ((1 1)) sum-%offset%))
  78           (setf (f2cl-lib:fref sum-%data% (k) ((1 1)) sum-%offset%) sigma)
  79           (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  80                         ((> i n) nil)
  81             (tagbody
  82               (setf sigma
  83                       (f2cl-lib:fref qr-%data%
  84                                      (i k)
  85                                      ((1 ndim) (1 1))
  86                                      qr-%offset%))
  87               (setf (f2cl-lib:fref qr-%data%
  88                                    (i k)
  89                                    ((1 ndim) (1 1))
  90                                    qr-%offset%)
  91                       (f2cl-lib:fref qr-%data%
  92                                      (i jbar)
  93                                      ((1 ndim) (1 1))
  94                                      qr-%offset%))
  95               (setf (f2cl-lib:fref qr-%data%
  96                                    (i jbar)
  97                                    ((1 ndim) (1 1))
  98                                    qr-%offset%)
  99                       sigma)
 100              label50))
 101          label70
 102           (setf sigma
 103                   (ddot (f2cl-lib:int-add (f2cl-lib:int-sub n k) 1)
 104                    (f2cl-lib:array-slice qr-%data%
 105                                          double-float
 106                                          (k k)
 107                                          ((1 ndim) (1 1))
 108                                          qr-%offset%)
 109                    1
 110                    (f2cl-lib:array-slice qr-%data%
 111                                          double-float
 112                                          (k k)
 113                                          ((1 ndim) (1 1))
 114                                          qr-%offset%)
 115                    1))
 116           (if (/= sigma 0.0f0) (go label60))
 117           (setf ierr 1)
 118           (go end_label)
 119          label60
 120           (if (= k n) (go label500))
 121           (setf qrkk
 122                   (f2cl-lib:fref qr-%data% (k k) ((1 ndim) (1 1)) qr-%offset%))
 123           (setf alphak (- (f2cl-lib:fsqrt sigma)))
 124           (if (< qrkk 0.0f0) (setf alphak (- alphak)))
 125           (setf (f2cl-lib:fref alpha-%data% (k) ((1 n)) alpha-%offset%) alphak)
 126           (setf beta (/ 1.0f0 (- sigma (* qrkk alphak))))
 127           (setf (f2cl-lib:fref qr-%data% (k k) ((1 ndim) (1 1)) qr-%offset%)
 128                   (- qrkk alphak))
 129           (f2cl-lib:fdo (j kp1 (f2cl-lib:int-add j 1))
 130                         ((> j np1) nil)
 131             (tagbody
 132              label80
 133               (setf (f2cl-lib:fref y-%data% (j) ((1 1)) y-%offset%)
 134                       (* beta
 135                          (ddot (f2cl-lib:int-add (f2cl-lib:int-sub n k) 1)
 136                           (f2cl-lib:array-slice qr-%data%
 137                                                 double-float
 138                                                 (k k)
 139                                                 ((1 ndim) (1 1))
 140                                                 qr-%offset%)
 141                           1
 142                           (f2cl-lib:array-slice qr-%data%
 143                                                 double-float
 144                                                 (k j)
 145                                                 ((1 ndim) (1 1))
 146                                                 qr-%offset%)
 147                           1)))))
 148           (f2cl-lib:fdo (j kp1 (f2cl-lib:int-add j 1))
 149                         ((> j np1) nil)
 150             (tagbody
 151               (f2cl-lib:fdo (i k (f2cl-lib:int-add i 1))
 152                             ((> i n) nil)
 153                 (tagbody
 154                   (setf (f2cl-lib:fref qr-%data%
 155                                        (i j)
 156                                        ((1 ndim) (1 1))
 157                                        qr-%offset%)
 158                           (-
 159                            (f2cl-lib:fref qr-%data%
 160                                           (i j)
 161                                           ((1 ndim) (1 1))
 162                                           qr-%offset%)
 163                            (*
 164                             (f2cl-lib:fref qr-%data%
 165                                            (i k)
 166                                            ((1 ndim) (1 1))
 167                                            qr-%offset%)
 168                             (f2cl-lib:fref y-%data% (j) ((1 1)) y-%offset%))))
 169                  label90))
 170               (setf (f2cl-lib:fref sum-%data% (j) ((1 1)) sum-%offset%)
 171                       (- (f2cl-lib:fref sum-%data% (j) ((1 1)) sum-%offset%)
 172                          (expt
 173                           (f2cl-lib:fref qr-%data%
 174                                          (k j)
 175                                          ((1 ndim) (1 1))
 176                                          qr-%offset%)
 177                           2)))
 178              label100))
 179          label500))
 180       (setf (f2cl-lib:fref alpha-%data% (n) ((1 n)) alpha-%offset%)
 181               (f2cl-lib:fref qr-%data% (n n) ((1 ndim) (1 1)) qr-%offset%))
 182       (go end_label)
 183      end_label
 184       (return (values nil nil nil nil nil ierr nil nil)))))
 185
 186 (in-package #:cl-user)
 187 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 188 (eval-when (:load-toplevel :compile-toplevel :execute)
 189   (setf (gethash 'fortran-to-lisp::dcpose
 190                  fortran-to-lisp::*f2cl-function-info*)
 191           (fortran-to-lisp::make-f2cl-finfo
 192            :arg-types '((fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 193                         (array double-float (*)) (array double-float (*))
 194                         (array fortran-to-lisp::integer4 (*))
 195                         (fortran-to-lisp::integer4) (array double-float (*))
 196                         (array double-float (*)))
 197            :return-values '(nil nil nil nil nil fortran-to-lisp::ierr nil nil)
 198            :calls '(fortran-to-lisp::ddot))))
 199