660 | | |
| 660 | |
| 661 | (defun print-zheev-results (e-val e-vec) |
| 662 | (format t "~2%ZHEEV Example Program Results~2%") |
| 663 | (let ((n (length e-val))) |
| 664 | (format t "Eigenvalues~%") |
| 665 | (dotimes (k n) |
| 666 | (format t " ~A" (aref e-val k))) |
| 667 | (format t "~2%Eigenvectors~%") |
| 668 | (dotimes (k n) |
| 669 | (format t "~%Eigenvector(~D)~%" k) |
| 670 | (dotimes (row n) |
| 671 | (format t "~A~%" (aref e-vec row k))) |
| 672 | (terpri)))) |
| 673 | |
| 674 | (defun test-zheev () |
| 675 | ;; The matrix is |
| 676 | ;; 1 #c(2 -1) #c(3 -1) #c(4 -1) |
| 677 | ;; #c(2 1) 2 #c(3 -2) #c(4 -2) |
| 678 | ;; #c(3 1) #c(3 2) 3 #c(4 -3) |
| 679 | ;; #c(4 1) #c(4 2) #c(4 3) 4 |
| 680 | ;; Recall that Fortran arrays are column-major order! |
| 681 | (let* ((n 4) |
| 682 | (a-mat (make-array (* n n) |
| 683 | :element-type '(complex double-float) |
| 684 | :initial-contents '(#c(1d0 0) |
| 685 | #c(2d0 1d0) |
| 686 | #c(3d0 1d0) |
| 687 | #c(4d0 1d0) |
| 688 | |
| 689 | #c(2d0 -1d0) |
| 690 | #c(2d0 0) |
| 691 | #c(3d0 2d0) |
| 692 | #c(4d0 2d0) |
| 693 | |
| 694 | #c(3d0 -1d0) |
| 695 | #c(3d0 -2d0) |
| 696 | #c(3d0 0) |
| 697 | #c(4d0 3d0) |
| 698 | |
| 699 | #c(4d0 -1d0) |
| 700 | #c(4d0 -2d0) |
| 701 | #c(4d0 -3d0) |
| 702 | #c(4d0 0)))) |
| 703 | (lwork 132) |
| 704 | (w (make-array n :element-type 'double-float)) |
| 705 | (work (make-array lwork :element-type '(complex double-float))) |
| 706 | (rwork (make-array (- (* 3 n) 2) :element-type 'double-float))) |
| 707 | (multiple-value-bind (z-jobz z-uplo z-n z-a z-lda z-w z-work z-lwork |
| 708 | z-rwork info) |
| 709 | (zheev "V" "U" n a-mat n w work lwork rwork 0) |
| 710 | (declare (ignore z-jobz z-uplo z-n z-a z-lda z-w z-work z-lwork |
| 711 | z-rwork)) |
| 712 | (cond ((zerop info) |
| 713 | (print-zheev-results w |
| 714 | (transpose (make-complex-eigvec n a-mat)))) |
| 715 | (t |
| 716 | (format t "Failure in ZHEEV. INFO = ~D~%" info))) |
| 717 | (format t "Optimum workspace required = ~D~%" (truncate (realpart (aref work 0)))) |
| 718 | (format t "Workspace provided = ~D~%" lwork) |
| 719 | (values w (transpose (make-complex-eigvec n a-mat)))))) |
| 720 | |
| 721 | (rt:deftest zheev.1 |
| 722 | (multiple-value-bind (val vec) |
| 723 | (test-zheev) |
| 724 | (let ((cval (make-array 4 :element-type '(complex double-float)))) |
| 725 | (map-into cval #'(lambda (x) |
| 726 | (complex x 0d0)) |
| 727 | val) |
| 728 | (list (check-eigen-val-vec 0 cval vec |
| 729 | #c(-4.244305402383179d0 0) |
| 730 | #(#C(0.38390086708134163d0 0.29405658581351474d0) |
| 731 | #C(0.4512081881756226d0 -0.11018120938621835d0) |
| 732 | #C(-0.026341826788003057d0 -0.48569831719444095d0) |
| 733 | #C(-0.5602011901552105d0 -0.0d0))) |
| 734 | (check-eigen-val-vec 1 cval vec |
| 735 | #c(-0.6885811461174348d0 0d0) |
| 736 | #(#C(-0.3975174811710425d0 0.5105011667992776d0) |
| 737 | #C(0.39532449113598805d0 -0.323828267501095d0) |
| 738 | #C(-0.43094682808194457d0 0.0382543464880426d0) |
| 739 | #C(0.36475148673610636d0 0.0d0))) |
| 740 | (check-eigen-val-vec 2 cval vec |
| 741 | #c(1.1412485214653287d0 0d0) |
| 742 | #(#C(0.37459176938191946d0 0.24136603676776178d0) |
| 743 | #C(-0.28951694726673083d0 0.4917390332315046d0) |
| 744 | #C(-0.37679299781495784d0 -0.3993997903501625d0) |
| 745 | #C(0.4174960446688179d0 0.0d0))) |
| 746 | (check-eigen-val-vec 3 cval vec |
| 747 | #c(13.791638027035287d0 0d0) |
| 748 | #(#C(0.3309009213351469d0 -0.19861339914465895d0) |
| 749 | #C(0.3727831818692677d0 -0.24193063397200412d0) |
| 750 | #C(0.4869962523459657d0 -0.19381997741123178d0) |
| 751 | #C(0.6154900747846208d0 0.0d0)))))) |
| 752 | (t t t t)) |
| 753 | |