Changeset 1155:d65e4d25ed88

Show
Ignore:
Timestamp:
11/13/13 19:23:07 (8 months ago)
Author:
Raymond Toy <toy.raymond@…>
Branch:
default
Tags:
tip
Message:

Add test for zheev.

Location:
packages
Files:
2 modified

Legend:

Unmodified
Added
Removed
  • packages/lapack/lapack-tests.lisp

    r1139 r1155  
    658658                                   #c(0.889824013759239d0 +0.000000000000000d0))))) 
    659659  (t t t t)) 
    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             
    661754(defun do-all-lapack-tests () 
    662755  (test-dgeev) 
     
    665758  (test-dgesdd) 
    666759  (test-dgesvd) 
    667   (test-zgeev)) 
     760  (test-zgeev) 
     761  (test-zheev)) 
    668762 
    669763;;; $Log$ 
  • packages/machar.lisp

    r1018 r1155  
    5050                the significand. 
    5151" 
     52  (declare (ignore ibeta it irnd ngrd machep negep iexp minexp  
     53                   maxexp eps epsneg xmin xmax)) 
    5254  (let ((ibeta (float-radix 1d0)) 
    5355        (it (float-digits 1d0)) 
     
    8890    (values ibeta it irnd ngrd machep negep iexp minexp  
    8991                        maxexp eps epsneg xmin xmax))) 
    90          
    91            
    92  
    93