home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SRC.DI$ / MEXTEST1.C < prev    next >
Encoding:
C/C++ Source or Header  |  1993-03-26  |  3.0 KB  |  162 lines
  1. #include <math.h>
  2. #include "mex.h"
  3.  
  4. /*
  5.  * MATLAB's MEX-file facility lets you call your own C routines
  6.  * directly from MATLAB.   Conversely, it is also possible to call
  7.  * internal MATLAB routines from your C routines, using the special
  8.  * 'mexCallMATLAB' utility routine.  In the way of an example, we've put
  9.  * together this routine, EIGS.C.
  10.  *
  11.  * This routine first forms and displays the matrix (in MATLAB notation):
  12.  *
  13.  *      hankel(1:4,4:-1:1) + sqrt(-1)*toeplitz(1:4,1:4)
  14.  *
  15.  * Next it finds the eigenvalues and eigenvectors (using the MATLAB
  16.  * function EIG), and displays the eigenvalue matrix.  Then it
  17.  * calculates the inverse of the eigenvalues to demonstrate manipulation
  18.  * of MATLAB results and how to deal with complex arithemetic.  Finally,
  19.  * the program displays the inverse values.
  20.  *
  21.  * Author: L. Shure 3-15-88
  22.  * Copyright 1988, The MathWorks, Inc., All rights reserved.
  23.  */
  24.  
  25. #define XR(i,j) xr[i+4*j]
  26. #define XI(i,j) xi[i+4*j]
  27.  
  28. static
  29. #ifdef __STDC__
  30. void fill_array(
  31.     double    *xr,
  32.     double    *xi
  33.     )
  34. #else
  35. fill_array(xr,xi)
  36. double *xr, *xi;
  37. #endif
  38. {
  39.     double tmp;
  40.     int i,j,jj;
  41. /*
  42.  * Remember, MATLAB stores matrices in their transposed form,
  43.  * i.e., columnwise, like FORTRAN.
  44.  *
  45.  * Fill real and imaginary parts of array.
  46.  */
  47.     for (j = 0; j < 4; j++) {
  48.         for (i = 0; i <= j; i++) {
  49.             XR(i,j) = 4 + i - j;
  50.             XR(j,i) = XR(i,j);
  51.             XI(i,j) = j - i + 1;
  52.             XI(j,i) = XI(i,j);
  53.         }
  54.     }
  55. /*
  56.  * Reorder columns of xr.
  57.  */
  58.     for (j = 0; j < 2; j++) {
  59.         for (i = 0; i < 4; i++) {
  60.             tmp = XR(i,j);
  61.             jj = 3 - j;
  62.             XR(i,j) = XR(i,jj);
  63.             XR(i,jj) = tmp;
  64.         }
  65.     }
  66. }
  67.  
  68. /*
  69.  * Invert diagonal elements of complex matrix of order 4
  70.  */
  71. static
  72. #ifdef __STDC__
  73. void invertd(
  74.     double    *xr,
  75.     double    *xi
  76.     )
  77. #else
  78. invertd(xr,xi)
  79. double *xr, *xi;
  80. #endif
  81. {
  82.     double tmp;
  83.     double *rx, *ix;
  84.     int i;
  85.  
  86.     rx = xr;
  87.     ix = xi;
  88. /*
  89.  * I know diagnonal elements of a 4 X 4 are 1:5:16
  90.  */
  91.     for (i = 0; i < 16; i += 5, rx += 5, ix += 5) {
  92.         tmp = *rx * *rx + *ix * *ix;
  93.         *rx = *rx / tmp;
  94.         *ix = - *ix / tmp;
  95.     }        
  96. }
  97. #ifdef __STDC__
  98. void
  99. mexFunction(
  100.     int        nlhs,
  101.     Matrix    *plhs[],
  102.     int        nrhs,
  103.     Matrix    *prhs[]
  104.     )
  105. #else
  106. mexFunction(nlhs, plhs, nrhs,prhs)
  107. int nlhs, nrhs;
  108. Matrix *plhs[], *prhs[];
  109. #endif
  110. {
  111.     int m, n, mrhs, mlhs;
  112.     Matrix *lhs[2], *x;
  113.  
  114.     m = n = 4;
  115.  
  116. /*
  117.  * allocate x matrix
  118.  */
  119.     x =  mxCreateFull(m, n, COMPLEX);
  120. /*
  121.  * create values in some arrays -- remember, MATLAB stores
  122.  * matrices column-wise
  123.  */
  124.     fill_array(mxGetPr(x), mxGetPi(x));
  125. /*
  126.  * print out initial matrix
  127.  */
  128.     mrhs = 1;
  129.     mlhs = 0;
  130.      if (!nlhs)
  131.         mexCallMATLAB(mlhs,lhs, mrhs, &x, "disp");
  132. /*
  133.  * calculate eigenvalues and eigenvectors
  134.  */
  135.     mlhs = 2;
  136.     mexCallMATLAB(mlhs, lhs, mrhs, &x, "eig");
  137. /*
  138.  * print out eigenvalue matrix
  139.  */
  140.     mlhs = 0;
  141.      if (!nlhs)
  142.         mexCallMATLAB(mlhs, lhs, mrhs, &lhs[1], "disp");
  143. /*
  144.  * take inverse of complex eigenvalues, just on diagonal
  145.  */
  146.     invertd(mxGetPr(lhs[1]), mxGetPi(lhs[1]));
  147. /*
  148.  * and print these out
  149.  */
  150.      if (!nlhs)
  151.         mexCallMATLAB(mlhs, lhs, mrhs, &lhs[1], "disp");
  152. /*
  153.  * Free allocated matrices
  154.  */
  155.     mxFreeMatrix(x);
  156.     mxFreeMatrix(lhs[1]);
  157.     if (nlhs)
  158.         plhs[0] = lhs[0];
  159.     else
  160.         mxFreeMatrix(lhs[0]);
  161. }
  162.