1;;; Compiled by f2cl version: 2;;; ("f2cl1.l,v 2edcbd958861 2012/05/30 03:34:52 toy $" 3;;; "f2cl2.l,v 96616d88fb7e 2008/02/22 22:19:34 rtoy $" 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 3fe93de3be82 2012/05/06 02:17:14 toy $" 7;;; "f2cl6.l,v 1d5cbacbb977 2008/08/24 00:56:27 rtoy $" 8;;; "macros.l,v 3fe93de3be82 2012/05/06 02:17:14 toy $") 9 10;;; Using Lisp CMU Common Lisp 20d (20D 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 :lapack) 18 19 20(let* ((zero 0.0) (one 1.0) (two 2.0) (const 1.5) (nsmax 15) (lds nsmax)) 21 (declare (type (double-float 0.0 0.0) zero) 22 (type (double-float 1.0 1.0) one) 23 (type (double-float 2.0 2.0) two) 24 (type (double-float 1.5 1.5) const) 25 (type (f2cl-lib:integer4 15 15) nsmax) 26 (type (f2cl-lib:integer4) lds) 27 (ignorable zero one two const nsmax lds)) 28 (defun dhseqr (job compz n ilo ihi h ldh wr wi z ldz work lwork info) 29 (declare (type (array double-float (*)) work z wi wr h) 30 (type (f2cl-lib:integer4) info lwork ldz ldh ihi ilo n) 31 (type (simple-string *) compz job)) 32 (f2cl-lib:with-multi-array-data 33 ((job character job-%data% job-%offset%) 34 (compz character compz-%data% compz-%offset%) 35 (h double-float h-%data% h-%offset%) 36 (wr double-float wr-%data% wr-%offset%) 37 (wi double-float wi-%data% wi-%offset%) 38 (z double-float z-%data% z-%offset%) 39 (work double-float work-%data% work-%offset%)) 40 (prog ((s 41 (make-array (the fixnum (reduce #'* (list lds nsmax))) 42 :element-type 'double-float)) 43 (v 44 (make-array (f2cl-lib:int-add nsmax 1) 45 :element-type 'double-float)) 46 (vv 47 (make-array (f2cl-lib:int-add nsmax 1) 48 :element-type 'double-float)) 49 (absw 0.0) (ovfl 0.0) (smlnum 0.0) (tau 0.0) (temp 0.0) (tst1 0.0) 50 (ulp 0.0) (unfl 0.0) (i 0) (i1 0) (i2 0) (ierr 0) (ii 0) (itemp 0) 51 (itn 0) (its 0) (j 0) (k 0) (l 0) (maxb 0) (nh 0) (nr 0) (ns 0) 52 (nv 0) (initz nil) (lquery nil) (wantt nil) (wantz nil)) 53 (declare (type (array double-float (*)) s v vv) 54 (type (double-float) absw ovfl smlnum tau temp tst1 ulp unfl) 55 (type (f2cl-lib:integer4) i i1 i2 ierr ii itemp itn its j k l 56 maxb nh nr ns nv) 57 (type f2cl-lib:logical initz lquery wantt wantz)) 58 (setf wantt (lsame job "S")) 59 (setf initz (lsame compz "I")) 60 (setf wantz (or initz (lsame compz "V"))) 61 (setf info 0) 62 (setf (f2cl-lib:fref work-%data% (1) ((1 *)) work-%offset%) 63 (coerce 64 (the f2cl-lib:integer4 65 (max (the f2cl-lib:integer4 1) 66 (the f2cl-lib:integer4 n))) 67 'double-float)) 68 (setf lquery (coerce (= lwork -1) 'f2cl-lib:logical)) 69 (cond 70 ((and (not (lsame job "E")) (not wantt)) 71 (setf info -1)) 72 ((and (not (lsame compz "N")) (not wantz)) 73 (setf info -2)) 74 ((< n 0) 75 (setf info -3)) 76 ((or (< ilo 1) 77 (> ilo 78 (max (the f2cl-lib:integer4 1) (the f2cl-lib:integer4 n)))) 79 (setf info -4)) 80 ((or 81 (< ihi (min (the f2cl-lib:integer4 ilo) (the f2cl-lib:integer4 n))) 82 (> ihi n)) 83 (setf info -5)) 84 ((< ldh (max (the f2cl-lib:integer4 1) (the f2cl-lib:integer4 n))) 85 (setf info -7)) 86 ((or (< ldz 1) 87 (and wantz 88 (< ldz 89 (max (the f2cl-lib:integer4 1) 90 (the f2cl-lib:integer4 n))))) 91 (setf info -11)) 92 ((and 93 (< lwork (max (the f2cl-lib:integer4 1) (the f2cl-lib:integer4 n))) 94 (not lquery)) 95 (setf info -13))) 96 (cond 97 ((/= info 0) 98 (xerbla "DHSEQR" (f2cl-lib:int-sub info)) 99 (go end_label)) 100 (lquery 101 (go end_label))) 102 (if initz (dlaset "Full" n n zero one z ldz)) 103 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 104 ((> i (f2cl-lib:int-add ilo (f2cl-lib:int-sub 1))) nil) 105 (tagbody 106 (setf (f2cl-lib:fref wr-%data% (i) ((1 *)) wr-%offset%) 107 (f2cl-lib:fref h-%data% (i i) ((1 ldh) (1 *)) h-%offset%)) 108 (setf (f2cl-lib:fref wi-%data% (i) ((1 *)) wi-%offset%) zero) 109 label10)) 110 (f2cl-lib:fdo (i (f2cl-lib:int-add ihi 1) (f2cl-lib:int-add i 1)) 111 ((> i n) nil) 112 (tagbody 113 (setf (f2cl-lib:fref wr-%data% (i) ((1 *)) wr-%offset%) 114 (f2cl-lib:fref h-%data% (i i) ((1 ldh) (1 *)) h-%offset%)) 115 (setf (f2cl-lib:fref wi-%data% (i) ((1 *)) wi-%offset%) zero) 116 label20)) 117 (if (= n 0) (go end_label)) 118 (cond 119 ((= ilo ihi) 120 (setf (f2cl-lib:fref wr-%data% (ilo) ((1 *)) wr-%offset%) 121 (f2cl-lib:fref h-%data% 122 (ilo ilo) 123 ((1 ldh) (1 *)) 124 h-%offset%)) 125 (setf (f2cl-lib:fref wi-%data% (ilo) ((1 *)) wi-%offset%) zero) 126 (go end_label))) 127 (f2cl-lib:fdo (j ilo (f2cl-lib:int-add j 1)) 128 ((> j (f2cl-lib:int-add ihi (f2cl-lib:int-sub 2))) nil) 129 (tagbody 130 (f2cl-lib:fdo (i (f2cl-lib:int-add j 2) (f2cl-lib:int-add i 1)) 131 ((> i n) nil) 132 (tagbody 133 (setf (f2cl-lib:fref h-%data% (i j) ((1 ldh) (1 *)) h-%offset%) 134 zero) 135 label30)) 136 label40)) 137 (setf nh (f2cl-lib:int-add (f2cl-lib:int-sub ihi ilo) 1)) 138 (setf ns (ilaenv 4 "DHSEQR" (f2cl-lib:f2cl-// job compz) n ilo ihi -1)) 139 (setf maxb 140 (ilaenv 8 "DHSEQR" (f2cl-lib:f2cl-// job compz) n ilo ihi -1)) 141 (cond 142 ((or (<= ns 2) (> ns nh) (>= maxb nh)) 143 (multiple-value-bind 144 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 145 var-10 var-11 var-12 var-13) 146 (dlahqr wantt wantz n ilo ihi h ldh wr wi ilo ihi z ldz info) 147 (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 148 var-8 var-9 var-10 var-11 var-12)) 149 (setf info var-13)) 150 (go end_label))) 151 (setf maxb 152 (max (the f2cl-lib:integer4 3) (the f2cl-lib:integer4 maxb))) 153 (setf ns 154 (min (the f2cl-lib:integer4 ns) 155 (the f2cl-lib:integer4 maxb) 156 (the f2cl-lib:integer4 nsmax))) 157 (setf unfl (dlamch "Safe minimum")) 158 (setf ovfl (/ one unfl)) 159 (multiple-value-bind (var-0 var-1) 160 (dlabad unfl ovfl) 161 (declare (ignore)) 162 (setf unfl var-0) 163 (setf ovfl var-1)) 164 (setf ulp (dlamch "Precision")) 165 (setf smlnum (* unfl (/ nh ulp))) 166 (cond 167 (wantt 168 (setf i1 1) 169 (setf i2 n))) 170 (setf itn (f2cl-lib:int-mul 30 nh)) 171 (setf i ihi) 172 label50 173 (setf l ilo) 174 (if (< i ilo) (go label170)) 175 (f2cl-lib:fdo (its 0 (f2cl-lib:int-add its 1)) 176 ((> its itn) nil) 177 (tagbody 178 (f2cl-lib:fdo (k i (f2cl-lib:int-add k (f2cl-lib:int-sub 1))) 179 ((> k (f2cl-lib:int-add l 1)) nil) 180 (tagbody 181 (setf tst1 182 (+ 183 (abs 184 (f2cl-lib:fref h-%data% 185 ((f2cl-lib:int-sub k 1) 186 (f2cl-lib:int-sub k 1)) 187 ((1 ldh) (1 *)) 188 h-%offset%)) 189 (abs 190 (f2cl-lib:fref h-%data% 191 (k k) 192 ((1 ldh) (1 *)) 193 h-%offset%)))) 194 (if (= tst1 zero) 195 (setf tst1 196 (dlanhs "1" 197 (f2cl-lib:int-add (f2cl-lib:int-sub i l) 1) 198 (f2cl-lib:array-slice h-%data% 199 double-float 200 (l l) 201 ((1 ldh) (1 *)) 202 h-%offset%) 203 ldh work))) 204 (if 205 (<= 206 (abs 207 (f2cl-lib:fref h-%data% 208 (k (f2cl-lib:int-sub k 1)) 209 ((1 ldh) (1 *)) 210 h-%offset%)) 211 (max (* ulp tst1) smlnum)) 212 (go label70)) 213 label60)) 214 label70 215 (setf l k) 216 (cond 217 ((> l ilo) 218 (setf (f2cl-lib:fref h-%data% 219 (l (f2cl-lib:int-sub l 1)) 220 ((1 ldh) (1 *)) 221 h-%offset%) 222 zero))) 223 (if (>= l (f2cl-lib:int-add (f2cl-lib:int-sub i maxb) 1)) 224 (go label160)) 225 (cond 226 ((not wantt) 227 (setf i1 l) 228 (setf i2 i))) 229 (cond 230 ((or (= its 20) (= its 30)) 231 (f2cl-lib:fdo (ii (f2cl-lib:int-add i (f2cl-lib:int-sub ns) 1) 232 (f2cl-lib:int-add ii 1)) 233 ((> ii i) nil) 234 (tagbody 235 (setf (f2cl-lib:fref wr-%data% (ii) ((1 *)) wr-%offset%) 236 (* const 237 (+ 238 (abs 239 (f2cl-lib:fref h-%data% 240 (ii (f2cl-lib:int-sub ii 1)) 241 ((1 ldh) (1 *)) 242 h-%offset%)) 243 (abs 244 (f2cl-lib:fref h-%data% 245 (ii ii) 246 ((1 ldh) (1 *)) 247 h-%offset%))))) 248 (setf (f2cl-lib:fref wi-%data% (ii) ((1 *)) wi-%offset%) 249 zero) 250 label80))) 251 (t 252 (dlacpy "Full" ns ns 253 (f2cl-lib:array-slice h-%data% 254 double-float 255 ((+ i (f2cl-lib:int-sub ns) 1) 256 (f2cl-lib:int-add 257 (f2cl-lib:int-sub i ns) 258 1)) 259 ((1 ldh) (1 *)) 260 h-%offset%) 261 ldh s lds) 262 (multiple-value-bind 263 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 264 var-9 var-10 var-11 var-12 var-13) 265 (dlahqr f2cl-lib:%false% f2cl-lib:%false% ns 1 ns s lds 266 (f2cl-lib:array-slice wr-%data% 267 double-float 268 ((+ i (f2cl-lib:int-sub ns) 1)) 269 ((1 *)) 270 wr-%offset%) 271 (f2cl-lib:array-slice wi-%data% 272 double-float 273 ((+ i (f2cl-lib:int-sub ns) 1)) 274 ((1 *)) 275 wi-%offset%) 276 1 ns z ldz ierr) 277 (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 278 var-7 var-8 var-9 var-10 var-11 var-12)) 279 (setf ierr var-13)) 280 (cond 281 ((> ierr 0) 282 (f2cl-lib:fdo (ii 1 (f2cl-lib:int-add ii 1)) 283 ((> ii ierr) nil) 284 (tagbody 285 (setf (f2cl-lib:fref wr-%data% 286 ((f2cl-lib:int-add 287 (f2cl-lib:int-sub i ns) 288 ii)) 289 ((1 *)) 290 wr-%offset%) 291 (f2cl-lib:fref s (ii ii) ((1 lds) (1 nsmax)))) 292 (setf (f2cl-lib:fref wi-%data% 293 ((f2cl-lib:int-add 294 (f2cl-lib:int-sub i ns) 295 ii)) 296 ((1 *)) 297 wi-%offset%) 298 zero) 299 label90)))))) 300 (setf (f2cl-lib:fref v (1) ((1 (f2cl-lib:int-add nsmax 1)))) one) 301 (f2cl-lib:fdo (ii 2 (f2cl-lib:int-add ii 1)) 302 ((> ii (f2cl-lib:int-add ns 1)) nil) 303 (tagbody 304 (setf (f2cl-lib:fref v (ii) ((1 (f2cl-lib:int-add nsmax 1)))) 305 zero) 306 label100)) 307 (setf nv 1) 308 (f2cl-lib:fdo (j (f2cl-lib:int-add i (f2cl-lib:int-sub ns) 1) 309 (f2cl-lib:int-add j 1)) 310 ((> j i) nil) 311 (tagbody 312 (cond 313 ((>= (f2cl-lib:fref wi (j) ((1 *))) zero) 314 (cond 315 ((= (f2cl-lib:fref wi (j) ((1 *))) zero) 316 (dcopy (f2cl-lib:int-add nv 1) v 1 vv 1) 317 (dgemv "No transpose" (f2cl-lib:int-add nv 1) nv one 318 (f2cl-lib:array-slice h-%data% 319 double-float 320 (l l) 321 ((1 ldh) (1 *)) 322 h-%offset%) 323 ldh vv 1 324 (- (f2cl-lib:fref wr-%data% (j) ((1 *)) wr-%offset%)) v 325 1) 326 (setf nv (f2cl-lib:int-add nv 1))) 327 ((> (f2cl-lib:fref wi (j) ((1 *))) zero) 328 (dcopy (f2cl-lib:int-add nv 1) v 1 vv 1) 329 (dgemv "No transpose" (f2cl-lib:int-add nv 1) nv one 330 (f2cl-lib:array-slice h-%data% 331 double-float 332 (l l) 333 ((1 ldh) (1 *)) 334 h-%offset%) 335 ldh v 1 336 (* (- two) 337 (f2cl-lib:fref wr-%data% (j) ((1 *)) wr-%offset%)) 338 vv 1) 339 (setf itemp (idamax (f2cl-lib:int-add nv 1) vv 1)) 340 (setf temp 341 (/ one 342 (max 343 (abs 344 (f2cl-lib:fref vv 345 (itemp) 346 ((1 347 (f2cl-lib:int-add nsmax 348 1))))) 349 smlnum))) 350 (dscal (f2cl-lib:int-add nv 1) temp vv 1) 351 (setf absw 352 (dlapy2 353 (f2cl-lib:fref wr-%data% 354 (j) 355 ((1 *)) 356 wr-%offset%) 357 (f2cl-lib:fref wi-%data% 358 (j) 359 ((1 *)) 360 wi-%offset%))) 361 (setf temp (* temp absw absw)) 362 (dgemv "No transpose" (f2cl-lib:int-add nv 2) 363 (f2cl-lib:int-add nv 1) one 364 (f2cl-lib:array-slice h-%data% 365 double-float 366 (l l) 367 ((1 ldh) (1 *)) 368 h-%offset%) 369 ldh vv 1 temp v 1) 370 (setf nv (f2cl-lib:int-add nv 2)))) 371 (setf itemp (idamax nv v 1)) 372 (setf temp 373 (abs 374 (f2cl-lib:fref v 375 (itemp) 376 ((1 (f2cl-lib:int-add nsmax 1)))))) 377 (cond 378 ((= temp zero) 379 (setf (f2cl-lib:fref v 380 (1) 381 ((1 (f2cl-lib:int-add nsmax 1)))) 382 one) 383 (f2cl-lib:fdo (ii 2 (f2cl-lib:int-add ii 1)) 384 ((> ii nv) nil) 385 (tagbody 386 (setf (f2cl-lib:fref v 387 (ii) 388 ((1 389 (f2cl-lib:int-add nsmax 1)))) 390 zero) 391 label110))) 392 (t 393 (setf temp (max temp smlnum)) 394 (dscal nv (/ one temp) v 1))))) 395 label120)) 396 (f2cl-lib:fdo (k l (f2cl-lib:int-add k 1)) 397 ((> k (f2cl-lib:int-add i (f2cl-lib:int-sub 1))) nil) 398 (tagbody 399 (setf nr 400 (min (the f2cl-lib:integer4 (f2cl-lib:int-add ns 1)) 401 (the f2cl-lib:integer4 402 (f2cl-lib:int-add (f2cl-lib:int-sub i k) 403 1)))) 404 (if (> k l) 405 (dcopy nr 406 (f2cl-lib:array-slice h-%data% 407 double-float 408 (k (f2cl-lib:int-sub k 1)) 409 ((1 ldh) (1 *)) 410 h-%offset%) 411 1 v 1)) 412 (multiple-value-bind (var-0 var-1 var-2 var-3 var-4) 413 (dlarfg nr 414 (f2cl-lib:fref v (1) ((1 (f2cl-lib:int-add nsmax 1)))) 415 (f2cl-lib:array-slice v 416 double-float 417 (2) 418 ((1 (f2cl-lib:int-add nsmax 1)))) 419 1 tau) 420 (declare (ignore var-0 var-2 var-3)) 421 (setf (f2cl-lib:fref v (1) ((1 (f2cl-lib:int-add nsmax 1)))) 422 var-1) 423 (setf tau var-4)) 424 (cond 425 ((> k l) 426 (setf (f2cl-lib:fref h-%data% 427 (k (f2cl-lib:int-sub k 1)) 428 ((1 ldh) (1 *)) 429 h-%offset%) 430 (f2cl-lib:fref v 431 (1) 432 ((1 (f2cl-lib:int-add nsmax 1))))) 433 (f2cl-lib:fdo (ii (f2cl-lib:int-add k 1) 434 (f2cl-lib:int-add ii 1)) 435 ((> ii i) nil) 436 (tagbody 437 (setf (f2cl-lib:fref h-%data% 438 (ii (f2cl-lib:int-sub k 1)) 439 ((1 ldh) (1 *)) 440 h-%offset%) 441 zero) 442 label130)))) 443 (setf (f2cl-lib:fref v (1) ((1 (f2cl-lib:int-add nsmax 1)))) 444 one) 445 (dlarfx "Left" nr (f2cl-lib:int-add (f2cl-lib:int-sub i2 k) 1) 446 v tau 447 (f2cl-lib:array-slice h-%data% 448 double-float 449 (k k) 450 ((1 ldh) (1 *)) 451 h-%offset%) 452 ldh work) 453 (dlarfx "Right" 454 (f2cl-lib:int-add 455 (f2cl-lib:int-sub 456 (min (the f2cl-lib:integer4 (f2cl-lib:int-add k nr)) 457 (the f2cl-lib:integer4 i)) 458 i1) 459 1) 460 nr v tau 461 (f2cl-lib:array-slice h-%data% 462 double-float 463 (i1 k) 464 ((1 ldh) (1 *)) 465 h-%offset%) 466 ldh work) 467 (cond 468 (wantz 469 (dlarfx "Right" nh nr v tau 470 (f2cl-lib:array-slice z-%data% 471 double-float 472 (ilo k) 473 ((1 ldz) (1 *)) 474 z-%offset%) 475 ldz work))) 476 label140)) 477 label150)) 478 (setf info i) 479 (go end_label) 480 label160 481 (multiple-value-bind 482 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 483 var-10 var-11 var-12 var-13) 484 (dlahqr wantt wantz n l i h ldh wr wi ilo ihi z ldz info) 485 (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 486 var-8 var-9 var-10 var-11 var-12)) 487 (setf info var-13)) 488 (if (> info 0) (go end_label)) 489 (setf itn (f2cl-lib:int-sub itn its)) 490 (setf i (f2cl-lib:int-sub l 1)) 491 (go label50) 492 label170 493 (setf (f2cl-lib:fref work-%data% (1) ((1 *)) work-%offset%) 494 (coerce 495 (the f2cl-lib:integer4 496 (max (the f2cl-lib:integer4 1) 497 (the f2cl-lib:integer4 n))) 498 'double-float)) 499 (go end_label) 500 end_label 501 (return 502 (values nil nil nil nil nil nil nil nil nil nil nil nil nil info)))))) 503 504(in-package #-gcl #:cl-user #+gcl "CL-USER") 505#+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or)) 506(eval-when (:load-toplevel :compile-toplevel :execute) 507 (setf (gethash 'fortran-to-lisp::dhseqr 508 fortran-to-lisp::*f2cl-function-info*) 509 (fortran-to-lisp::make-f2cl-finfo 510 :arg-types '((simple-string) (simple-string) 511 (fortran-to-lisp::integer4) (fortran-to-lisp::integer4) 512 (fortran-to-lisp::integer4) (array double-float (*)) 513 (fortran-to-lisp::integer4) (array double-float (*)) 514 (array double-float (*)) (array double-float (*)) 515 (fortran-to-lisp::integer4) (array double-float (*)) 516 (fortran-to-lisp::integer4) 517 (fortran-to-lisp::integer4)) 518 :return-values '(nil nil nil nil nil nil nil nil nil nil nil nil nil 519 fortran-to-lisp::info) 520 :calls '(fortran-to-lisp::dlarfx fortran-to-lisp::dlarfg 521 fortran-to-lisp::dlapy2 fortran-to-lisp::dscal 522 fortran-to-lisp::idamax fortran-to-lisp::dgemv 523 fortran-to-lisp::dcopy fortran-to-lisp::dlacpy 524 fortran-to-lisp::dlanhs fortran-to-lisp::dlabad 525 fortran-to-lisp::dlamch fortran-to-lisp::dlahqr 526 fortran-to-lisp::ilaenv fortran-to-lisp::dlaset 527 fortran-to-lisp::xerbla fortran-to-lisp::lsame)))) 528 529