Changeset 1155:d65e4d25ed88


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

Add test for zheev.

Location:
packages
Files:
2 edited

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  
Note: See TracChangeset for help on using the changeset viewer.