home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SRC.DI$ / MEXTST1.F < prev    next >
Encoding:
Text File  |  1993-03-26  |  4.1 KB  |  141 lines
  1.       subroutine EIGS(nlhs, plhs)
  2.       integer nlhs
  3.       integer plhs(*)
  4.  
  5. c MATLAB's MEX-file facility lets you call your own Fortran subroutines
  6. c directly from MATLAB.   Conversely, it is also possible to call
  7. c internal MATLAB routines from your Fortran routines, using the special
  8. c mexCallMATLAB utility routine.  In the way of an example, we've put together 
  9. c this subroutine, MEXTST1.F, together with its gateway routine, MEXTST1G.F.
  10.  
  11. c This routine first forms and displays the matrix (in MATLAB notation):
  12. c
  13. c      hankel(1:4,4:-1:1) + sqrt(-1)*toeplitz(1:4,1:4)
  14. c
  15. c Next it finds the eigenvalues and eigenvectors (using the
  16. c MATLAB function EIG), and displays the eigenvalue matrix.  Then it
  17. c calculates the inverse of the eigenvalues to demonstrate manipulation
  18. c of MATLAB results and how to deal with complex arithemetic.  Finally,
  19. c the program displays the inverse values.
  20.  
  21. c Author: L. Shure 3-15-88
  22. c Copyright 1988, The MathWorks, Inc., All rights reserved.
  23.  
  24. c Dimension arrays
  25.       double precision xreal(4,4), ximag(4,4)
  26.  
  27. c Declare MEX utility functions (see MATLAB user's guide)
  28.       integer mxCreateFull, mxGetPr, mxGetPi
  29.  
  30. c Declare all the variables we intend to use for the example.
  31.       integer m, n, imagflag, lhs(2), mlhs, mrhs, x, xr, xi, mn, rhs
  32.       integer evaluer, evaluei
  33.  
  34. c Fill xreal and ximag arrays with the matrix (in MATLAB notation)
  35. c hankel(1:4,4:-1:1)+sqrt(-1)*toeplitz(1:4,1:4)
  36.       call fill(xreal,ximag)
  37.       m=4
  38.       n=4
  39.       mn = m*n
  40.       imagflag = 1
  41.  
  42. c Create (allocate) a complex matrix of size m-by-n
  43.       x = MXCREATEFULL(m,n,imagflag)
  44.  
  45. c Get pointers to the real and imaginary parts of the matrix
  46.       xr = MXGETPR(x)
  47.       xi = MXGETPI(x)
  48.  
  49. c Copy our arrays xreal and ximag to their corresponding places in
  50. c the matrix structure x.
  51.       call MXCOPYREAL8TOPTR(xreal,xr,mn)
  52.       call MXCOPYREAL8TOPTR(ximag,xi,mn)
  53.  
  54. c Find eigenvectors and eigenvalues of matrix x using the MATLAB
  55. c function EIG.  We give mexCallMATLAB one right-hand-argument and expect
  56. c two left-hand-side arrays to be returned.
  57.       mrhs = 1
  58.       mlhs = 2
  59.       call MEXCALLMATLAB(mlhs,lhs,mrhs,x,'eig')
  60.  
  61. c Display original x matrix, therefore there is no left-hand side.
  62.       mlhs = 0
  63.       if (nlhs .eq. 0) call MEXCALLMATLAB(mlhs,lhs,mrhs,x,'disp')
  64.  
  65. c Display eigenvalue matrix.
  66.       rhs = lhs(2)
  67.       if (nlhs .eq. 0) call MEXCALLMATLAB(mlhs,lhs,mrhs,rhs,'disp')
  68.  
  69. c Invert diagonal elements of the complex matrix by two methods.  The
  70. c first uses auxiliary arrays for performing the inversion, the second
  71. c does the calculations in place.
  72.       evaluer = MXGETPR(lhs(2))
  73.       evaluei = MXGETPI(lhs(2))
  74.  
  75. c Copy these to new working arrays.
  76.       call MXCOPYPTRTOREAL8(evaluer, xreal, mn)
  77.       call MXCOPYPTRTOREAL8(evaluei, ximag, mn)
  78.       call invertd(xreal, ximag)
  79.  
  80. c Copy answers back to a matrix array x
  81.       call MXCOPYREAL8TOPTR(xreal,xr,mn)
  82.       call MXCOPYREAL8TOPTR(ximag,xi,mn)
  83.       if (nlhs .eq. 0) call MEXCALLMATLAB(mlhs,lhs,mrhs,x,'disp')
  84.  
  85. c Free all allocated matrices
  86.       call MXFREEMATRIX(x)
  87.       call MXFREEMATRIX(lhs(2))
  88.       if (nlhs .ne. 0) goto 1000
  89.       call MXFREEMATRIX(lhs(1))
  90.       goto 1010
  91.  1000 plhs(1) = lhs(1)
  92.  1010 return
  93.       end
  94.  
  95.       subroutine fill(xr,xi)
  96.       double precision xr(4,4) ,xi(4,4) ,tmp
  97.       integer i ,j
  98.  
  99. c Form the matrix: hankel(1:4,4:-1:1)+sqrt(-1)*toeplitz(1:4,1:4)
  100.       do 1000 j=1,4
  101.       do 1000 i=1,j
  102.         xr(i,j) = 4+i-j
  103.         xr(j,i) = xr(i,j)
  104.         xi(i,j) = j-i+1
  105.  1000   xi(j,i) = xi(i,j)
  106.  
  107. c Reorder columns of xr.
  108.       do 1050 j=1,2
  109.       do 1050 i=1,4
  110.         tmp = xr(i,j)
  111.         jj = 5-j
  112.         xr(i,j) = xr(i,jj)
  113.  1050   xr(i,jj) = tmp
  114.       return
  115.       end
  116.                                          
  117.       subroutine invertd(xr,xi)
  118.       double precision xr(*), xi(*), tmp
  119.       integer i
  120. c
  121. c Just invert diagonal elements of a complex matrix of order 4.
  122. c
  123.       do 1000 i=1,16,5
  124.         tmp = xr(i)**2+xi(i)**2
  125.         xr(i) = xr(i)/tmp
  126.  1000   xi(i) = -xi(i)/tmp
  127.       return
  128.       end
  129.  
  130.       subroutine copyrl(x,y,n)
  131. c
  132. c Copy one length n real array to another
  133. c
  134.       double precision x(*), y(*)
  135.       integer n, i
  136.                              
  137.       do 1000 i=1,n
  138.  1000   y(i) = x(i)
  139.       return
  140.       end
  141.