home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 1 / RISC_DISC_1.iso / pd_share / code / meschach / !Meschach / c / zmatop < prev    next >
Encoding:
Text File  |  1994-01-12  |  14.9 KB  |  613 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  5. **
  6. **                 Meschach Library
  7. ** 
  8. ** This Meschach Library is provided "as is" without any express 
  9. ** or implied warranty of any kind with respect to this software. 
  10. ** In particular the authors shall not be liable for any direct, 
  11. ** indirect, special, incidental or consequential damages arising 
  12. ** in any way from use of the software.
  13. ** 
  14. ** Everyone is granted permission to copy, modify and redistribute this
  15. ** Meschach Library, provided:
  16. **  1.  All copies contain this copyright notice.
  17. **  2.  All modified copies shall carry a notice stating who
  18. **      made the last modification and the date of such modification.
  19. **  3.  No charge is made for this software or works derived from it.  
  20. **      This clause shall not be construed as constraining other software
  21. **      distributed on the same medium as this software, nor is a
  22. **      distribution fee considered a charge.
  23. **
  24. ***************************************************************************/
  25.  
  26.  
  27.  
  28. #include    <stdio.h>
  29. #include    "zmatrix.h"
  30.  
  31. static    char    rcsid[] = "$Id: zmatop.c,v 1.1 1994/01/13 04:24:46 des Exp $";
  32.  
  33.  
  34. #define    is_zero(z)    ((z).re == 0.0 && (z).im == 0.0)
  35.  
  36. /* zm_add -- matrix addition -- may be in-situ */
  37. ZMAT    *zm_add(mat1,mat2,out)
  38. ZMAT    *mat1,*mat2,*out;
  39. {
  40.     u_int    m,n,i;
  41.     
  42.     if ( mat1==ZMNULL || mat2==ZMNULL )
  43.     error(E_NULL,"zm_add");
  44.     if ( mat1->m != mat2->m || mat1->n != mat2->n )
  45.     error(E_SIZES,"zm_add");
  46.     if ( out==ZMNULL || out->m != mat1->m || out->n != mat1->n )
  47.     out = zm_resize(out,mat1->m,mat1->n);
  48.     m = mat1->m;    n = mat1->n;
  49.     for ( i=0; i<m; i++ )
  50.     {
  51.     __zadd__(mat1->me[i],mat2->me[i],out->me[i],(int)n);
  52.     /**************************************************
  53.       for ( j=0; j<n; j++ )
  54.       out->me[i][j] = mat1->me[i][j]+mat2->me[i][j];
  55.       **************************************************/
  56.     }
  57.     
  58.     return (out);
  59. }
  60.  
  61. /* zm_sub -- matrix subtraction -- may be in-situ */
  62. ZMAT    *zm_sub(mat1,mat2,out)
  63. ZMAT    *mat1,*mat2,*out;
  64. {
  65.     u_int    m,n,i;
  66.     
  67.     if ( mat1==ZMNULL || mat2==ZMNULL )
  68.     error(E_NULL,"zm_sub");
  69.     if ( mat1->m != mat2->m || mat1->n != mat2->n )
  70.     error(E_SIZES,"zm_sub");
  71.     if ( out==ZMNULL || out->m != mat1->m || out->n != mat1->n )
  72.     out = zm_resize(out,mat1->m,mat1->n);
  73.     m = mat1->m;    n = mat1->n;
  74.     for ( i=0; i<m; i++ )
  75.     {
  76.     __zsub__(mat1->me[i],mat2->me[i],out->me[i],(int)n);
  77.     /**************************************************
  78.       for ( j=0; j<n; j++ )
  79.       out->me[i][j] = mat1->me[i][j]-mat2->me[i][j];
  80.     **************************************************/
  81.     }
  82.     
  83.     return (out);
  84. }
  85.  
  86. /*
  87.   Note: In the following routines, "adjoint" means complex conjugate
  88.   transpose:
  89.   A* = conjugate(A^T)
  90.   */
  91.  
  92. /* zm_mlt -- matrix-matrix multiplication */
  93. ZMAT    *zm_mlt(A,B,OUT)
  94. ZMAT    *A,*B,*OUT;
  95. {
  96.     u_int    i, /* j, */ k, m, n, p;
  97.     complex    **A_v, **B_v /*, *B_row, *OUT_row, sum, tmp */;
  98.     
  99.     if ( A==ZMNULL || B==ZMNULL )
  100.     error(E_NULL,"zm_mlt");
  101.     if ( A->n != B->m )
  102.     error(E_SIZES,"zm_mlt");
  103.     if ( A == OUT || B == OUT )
  104.     error(E_INSITU,"zm_mlt");
  105.     m = A->m;    n = A->n;    p = B->n;
  106.     A_v = A->me;        B_v = B->me;
  107.     
  108.     if ( OUT==ZMNULL || OUT->m != A->m || OUT->n != B->n )
  109.     OUT = zm_resize(OUT,A->m,B->n);
  110.     
  111.     /****************************************************************
  112.       for ( i=0; i<m; i++ )
  113.       for  ( j=0; j<p; j++ )
  114.       {
  115.       sum = 0.0;
  116.       for ( k=0; k<n; k++ )
  117.       sum += A_v[i][k]*B_v[k][j];
  118.       OUT->me[i][j] = sum;
  119.       }
  120.     ****************************************************************/
  121.     zm_zero(OUT);
  122.     for ( i=0; i<m; i++ )
  123.     for ( k=0; k<n; k++ )
  124.     {
  125.         if ( ! is_zero(A_v[i][k]) )
  126.         __zmltadd__(OUT->me[i],B_v[k],A_v[i][k],(int)p,Z_NOCONJ);
  127.         /**************************************************
  128.           B_row = B_v[k];    OUT_row = OUT->me[i];
  129.           for ( j=0; j<p; j++ )
  130.           (*OUT_row++) += tmp*(*B_row++);
  131.         **************************************************/
  132.     }
  133.     
  134.     return OUT;
  135. }
  136.  
  137. /* zmma_mlt -- matrix-matrix adjoint multiplication
  138.    -- A.B* is returned, and stored in OUT */
  139. ZMAT    *zmma_mlt(A,B,OUT)
  140. ZMAT    *A, *B, *OUT;
  141. {
  142.     int    i, j, limit;
  143.     /* complex    *A_row, *B_row, sum; */
  144.     
  145.     if ( ! A || ! B )
  146.     error(E_NULL,"zmma_mlt");
  147.     if ( A == OUT || B == OUT )
  148.     error(E_INSITU,"zmma_mlt");
  149.     if ( A->n != B->n )
  150.     error(E_SIZES,"zmma_mlt");
  151.     if ( ! OUT || OUT->m != A->m || OUT->n != B->m )
  152.     OUT = zm_resize(OUT,A->m,B->m);
  153.     
  154.     limit = A->n;
  155.     for ( i = 0; i < A->m; i++ )
  156.     for ( j = 0; j < B->m; j++ )
  157.     {
  158.         OUT->me[i][j] = __zip__(B->me[j],A->me[i],(int)limit,Z_CONJ);
  159.         /**************************************************
  160.           sum = 0.0;
  161.           A_row = A->me[i];
  162.           B_row = B->me[j];
  163.           for ( k = 0; k < limit; k++ )
  164.           sum += (*A_row++)*(*B_row++);
  165.           OUT->me[i][j] = sum;
  166.           **************************************************/
  167.     }
  168.     
  169.     return OUT;
  170. }
  171.  
  172. /* zmam_mlt -- matrix adjoint-matrix multiplication
  173.    -- A*.B is returned, result stored in OUT */
  174. ZMAT    *zmam_mlt(A,B,OUT)
  175. ZMAT    *A, *B, *OUT;
  176. {
  177.     int    i, k, limit;
  178.     /* complex    *B_row, *OUT_row, multiplier; */
  179.     complex    tmp;
  180.     
  181.     if ( ! A || ! B )
  182.     error(E_NULL,"zmam_mlt");
  183.     if ( A == OUT || B == OUT )
  184.     error(E_INSITU,"zmam_mlt");
  185.     if ( A->m != B->m )
  186.     error(E_SIZES,"zmam_mlt");
  187.     if ( ! OUT || OUT->m != A->n || OUT->n != B->n )
  188.     OUT = zm_resize(OUT,A->n,B->n);
  189.     
  190.     limit = B->n;
  191.     zm_zero(OUT);
  192.     for ( k = 0; k < A->m; k++ )
  193.     for ( i = 0; i < A->n; i++ )
  194.     {
  195.         tmp.re =   A->me[k][i].re;
  196.         tmp.im = - A->me[k][i].im;
  197.         if ( ! is_zero(tmp) )
  198.         __zmltadd__(OUT->me[i],B->me[k],tmp,(int)limit,Z_NOCONJ);
  199.     }
  200.     
  201.     return OUT;
  202. }
  203.  
  204. /* zmv_mlt -- matrix-vector multiplication 
  205.    -- Note: b is treated as a column vector */
  206. ZVEC    *zmv_mlt(A,b,out)
  207. ZMAT    *A;
  208. ZVEC    *b,*out;
  209. {
  210.     u_int    i, m, n;
  211.     complex    **A_v, *b_v /*, *A_row */;
  212.     /* register complex    sum; */
  213.     
  214.     if ( A==ZMNULL || b==ZVNULL )
  215.     error(E_NULL,"zmv_mlt");
  216.     if ( A->n != b->dim )
  217.     error(E_SIZES,"zmv_mlt");
  218.     if ( b == out )
  219.     error(E_INSITU,"zmv_mlt");
  220.     if ( out == ZVNULL || out->dim != A->m )
  221.     out = zv_resize(out,A->m);
  222.     
  223.     m = A->m;        n = A->n;
  224.     A_v = A->me;        b_v = b->ve;
  225.     for ( i=0; i<m; i++ )
  226.     {
  227.     /* for ( j=0; j<n; j++ )
  228.        sum += A_v[i][j]*b_v[j]; */
  229.     out->ve[i] = __zip__(A_v[i],b_v,(int)n,Z_NOCONJ);
  230.     /**************************************************
  231.       A_row = A_v[i];        b_v = b->ve;
  232.       for ( j=0; j<n; j++ )
  233.       sum += (*A_row++)*(*b_v++);
  234.       out->ve[i] = sum;
  235.     **************************************************/
  236.     }
  237.     
  238.     return out;
  239. }
  240.  
  241. /* zsm_mlt -- scalar-matrix multiply -- may be in-situ */
  242. ZMAT    *zsm_mlt(scalar,matrix,out)
  243. complex    scalar;
  244. ZMAT    *matrix,*out;
  245. {
  246.     u_int    m,n,i;
  247.     
  248.     if ( matrix==ZMNULL )
  249.     error(E_NULL,"zsm_mlt");
  250.     if ( out==ZMNULL || out->m != matrix->m || out->n != matrix->n )
  251.     out = zm_resize(out,matrix->m,matrix->n);
  252.     m = matrix->m;    n = matrix->n;
  253.     for ( i=0; i<m; i++ )
  254.     __zmlt__(matrix->me[i],scalar,out->me[i],(int)n);
  255.     /**************************************************
  256.       for ( j=0; j<n; j++ )
  257.       out->me[i][j] = scalar*matrix->me[i][j];
  258.       **************************************************/
  259.     return (out);
  260. }
  261.  
  262. /* zvm_mlt -- vector adjoint-matrix multiplication */
  263. ZVEC    *zvm_mlt(A,b,out)
  264. ZMAT    *A;
  265. ZVEC    *b,*out;
  266. {
  267.     u_int    j,m,n;
  268.     /* complex    sum,**A_v,*b_v; */
  269.     
  270.     if ( A==ZMNULL || b==ZVNULL )
  271.     error(E_NULL,"zvm_mlt");
  272.     if ( A->m != b->dim )
  273.     error(E_SIZES,"zvm_mlt");
  274.     if ( b == out )
  275.     error(E_INSITU,"zvm_mlt");
  276.     if ( out == ZVNULL || out->dim != A->n )
  277.     out = zv_resize(out,A->n);
  278.     
  279.     m = A->m;        n = A->n;
  280.     
  281.     zv_zero(out);
  282.     for ( j = 0; j < m; j++ )
  283.     if ( b->ve[j].re != 0.0 || b->ve[j].im != 0.0  )
  284.         __zmltadd__(out->ve,A->me[j],b->ve[j],(int)n,Z_CONJ);
  285.     /**************************************************
  286.       A_v = A->me;        b_v = b->ve;
  287.       for ( j=0; j<n; j++ )
  288.       {
  289.       sum = 0.0;
  290.       for ( i=0; i<m; i++ )
  291.       sum += b_v[i]*A_v[i][j];
  292.       out->ve[j] = sum;
  293.       }
  294.       **************************************************/
  295.     
  296.     return out;
  297. }
  298.  
  299. /* zm_adjoint -- adjoint matrix */
  300. ZMAT    *zm_adjoint(in,out)
  301. ZMAT    *in, *out;
  302. {
  303.     int    i, j;
  304.     int    in_situ;
  305.     complex    tmp;
  306.     
  307.     if ( in == ZMNULL )
  308.     error(E_NULL,"zm_adjoint");
  309.     if ( in == out && in->n != in->m )
  310.     error(E_INSITU2,"zm_adjoint");
  311.     in_situ = ( in == out );
  312.     if ( out == ZMNULL || out->m != in->n || out->n != in->m )
  313.     out = zm_resize(out,in->n,in->m);
  314.     
  315.     if ( ! in_situ )
  316.     {
  317.     for ( i = 0; i < in->m; i++ )
  318.         for ( j = 0; j < in->n; j++ )
  319.         {
  320.         out->me[j][i].re =   in->me[i][j].re;
  321.         out->me[j][i].im = - in->me[i][j].im;
  322.         }
  323.     }
  324.     else
  325.     {
  326.     for ( i = 0 ; i < in->m; i++ )
  327.     {
  328.         for ( j = 0; j < i; j++ )
  329.         {
  330.         tmp.re = in->me[i][j].re;
  331.         tmp.im = in->me[i][j].im;
  332.         in->me[i][j].re =   in->me[j][i].re;
  333.         in->me[i][j].im = - in->me[j][i].im;
  334.         in->me[j][i].re =   tmp.re;
  335.         in->me[j][i].im = - tmp.im;
  336.         }
  337.         in->me[i][i].im = - in->me[i][i].im;
  338.     }
  339.     }
  340.     
  341.     return out;
  342. }
  343.  
  344. /* zswap_rows -- swaps rows i and j of matrix A upto column lim */
  345. ZMAT    *zswap_rows(A,i,j,lo,hi)
  346. ZMAT    *A;
  347. int    i, j, lo, hi;
  348. {
  349.     int    k;
  350.     complex    **A_me, tmp;
  351.     
  352.     if ( ! A )
  353.     error(E_NULL,"swap_rows");
  354.     if ( i < 0 || j < 0 || i >= A->m || j >= A->m )
  355.     error(E_SIZES,"swap_rows");
  356.     lo = max(0,lo);
  357.     hi = min(hi,A->n-1);
  358.     A_me = A->me;
  359.     
  360.     for ( k = lo; k <= hi; k++ )
  361.     {
  362.     tmp = A_me[k][i];
  363.     A_me[k][i] = A_me[k][j];
  364.     A_me[k][j] = tmp;
  365.     }
  366.     return A;
  367. }
  368.  
  369. /* zswap_cols -- swap columns i and j of matrix A upto row lim */
  370. ZMAT    *zswap_cols(A,i,j,lo,hi)
  371. ZMAT    *A;
  372. int    i, j, lo, hi;
  373. {
  374.     int    k;
  375.     complex    **A_me, tmp;
  376.     
  377.     if ( ! A )
  378.     error(E_NULL,"swap_cols");
  379.     if ( i < 0 || j < 0 || i >= A->n || j >= A->n )
  380.     error(E_SIZES,"swap_cols");
  381.     lo = max(0,lo);
  382.     hi = min(hi,A->m-1);
  383.     A_me = A->me;
  384.     
  385.     for ( k = lo; k <= hi; k++ )
  386.     {
  387.     tmp = A_me[i][k];
  388.     A_me[i][k] = A_me[j][k];
  389.     A_me[j][k] = tmp;
  390.     }
  391.     return A;
  392. }
  393.  
  394. /* mz_mltadd -- matrix-scalar multiply and add
  395.    -- may be in situ
  396.    -- returns out == A1 + s*A2 */
  397. ZMAT    *mz_mltadd(A1,A2,s,out)
  398. ZMAT    *A1, *A2, *out;
  399. complex    s;
  400. {
  401.     /* register complex    *A1_e, *A2_e, *out_e; */
  402.     /* register int    j; */
  403.     int    i, m, n;
  404.     
  405.     if ( ! A1 || ! A2 )
  406.     error(E_NULL,"mz_mltadd");
  407.     if ( A1->m != A2->m || A1->n != A2->n )
  408.     error(E_SIZES,"mz_mltadd");
  409.     
  410.     if ( s.re == 0.0 && s.im == 0.0 )
  411.     return zm_copy(A1,out);
  412.     if ( s.re == 1.0 && s.im == 0.0 )
  413.     return zm_add(A1,A2,out);
  414.     
  415.     tracecatch(out = zm_copy(A1,out),"mz_mltadd");
  416.     
  417.     m = A1->m;    n = A1->n;
  418.     for ( i = 0; i < m; i++ )
  419.     {
  420.     __zmltadd__(out->me[i],A2->me[i],s,(int)n,Z_NOCONJ);
  421.     /**************************************************
  422.       A1_e = A1->me[i];
  423.       A2_e = A2->me[i];
  424.       out_e = out->me[i];
  425.       for ( j = 0; j < n; j++ )
  426.       out_e[j] = A1_e[j] + s*A2_e[j];
  427.       **************************************************/
  428.     }
  429.     
  430.     return out;
  431. }
  432.  
  433. /* zmv_mltadd -- matrix-vector multiply and add
  434.    -- may not be in situ
  435.    -- returns out == v1 + alpha*A*v2 */
  436. ZVEC    *zmv_mltadd(v1,v2,A,alpha,out)
  437. ZVEC    *v1, *v2, *out;
  438. ZMAT    *A;
  439. complex    alpha;
  440. {
  441.     /* register    int    j; */
  442.     int    i, m, n;
  443.     complex    tmp, *v2_ve, *out_ve;
  444.     
  445.     if ( ! v1 || ! v2 || ! A )
  446.     error(E_NULL,"zmv_mltadd");
  447.     if ( out == v2 )
  448.     error(E_INSITU,"zmv_mltadd");
  449.     if ( v1->dim != A->m || v2->dim != A-> n )
  450.     error(E_SIZES,"zmv_mltadd");
  451.     
  452.     tracecatch(out = zv_copy(v1,out),"zmv_mltadd");
  453.     
  454.     v2_ve = v2->ve;    out_ve = out->ve;
  455.     m = A->m;    n = A->n;
  456.     
  457.     if ( alpha.re == 0.0 && alpha.im == 0.0 )
  458.     return out;
  459.     
  460.     for ( i = 0; i < m; i++ )
  461.     {
  462.     tmp = __zip__(A->me[i],v2_ve,(int)n,Z_NOCONJ);
  463.     out_ve[i].re += alpha.re*tmp.re - alpha.im*tmp.im;
  464.     out_ve[i].im += alpha.re*tmp.im + alpha.im*tmp.re;
  465.     /**************************************************
  466.       A_e = A->me[i];
  467.       sum = 0.0;
  468.       for ( j = 0; j < n; j++ )
  469.       sum += A_e[j]*v2_ve[j];
  470.       out_ve[i] = v1->ve[i] + alpha*sum;
  471.       **************************************************/
  472.     }
  473.     
  474.     return out;
  475. }
  476.  
  477. /* zvm_mltadd -- vector-matrix multiply and add a la zvm_mlt()
  478.    -- may not be in situ
  479.    -- returns out == v1 + v2*.A */
  480. ZVEC    *zvm_mltadd(v1,v2,A,alpha,out)
  481. ZVEC    *v1, *v2, *out;
  482. ZMAT    *A;
  483. complex    alpha;
  484. {
  485.     int    /* i, */ j, m, n;
  486.     complex    tmp, /* *A_e, */ *out_ve;
  487.     
  488.     if ( ! v1 || ! v2 || ! A )
  489.     error(E_NULL,"zvm_mltadd");
  490.     if ( v2 == out )
  491.     error(E_INSITU,"zvm_mltadd");
  492.     if ( v1->dim != A->n || A->m != v2->dim )
  493.     error(E_SIZES,"zvm_mltadd");
  494.     
  495.     tracecatch(out = zv_copy(v1,out),"zvm_mltadd");
  496.     
  497.     out_ve = out->ve;    m = A->m;    n = A->n;
  498.     for ( j = 0; j < m; j++ )
  499.     {
  500.     /* tmp = zmlt(v2->ve[j],alpha); */
  501.     tmp.re =   v2->ve[j].re*alpha.re - v2->ve[j].im*alpha.im;
  502.     tmp.im =   v2->ve[j].re*alpha.im + v2->ve[j].im*alpha.re;
  503.     if ( tmp.re != 0.0 || tmp.im != 0.0 )
  504.         __zmltadd__(out_ve,A->me[j],tmp,(int)n,Z_CONJ);
  505.     /**************************************************
  506.       A_e = A->me[j];
  507.       for ( i = 0; i < n; i++ )
  508.       out_ve[i] += A_e[i]*tmp;
  509.     **************************************************/
  510.     }
  511.     
  512.     return out;
  513. }
  514.  
  515. /* zget_col -- gets a specified column of a matrix; returned as a vector */
  516. ZVEC    *zget_col(mat,col,vec)
  517. int    col;
  518. ZMAT    *mat;
  519. ZVEC    *vec;
  520. {
  521.     u_int    i;
  522.  
  523.     if ( mat==ZMNULL )
  524.         error(E_NULL,"zget_col");
  525.     if ( col < 0 || col >= mat->n )
  526.         error(E_RANGE,"zget_col");
  527.     if ( vec==ZVNULL || vec->dim<mat->m )
  528.         vec = zv_resize(vec,mat->m);
  529.  
  530.     for ( i=0; i<mat->m; i++ )
  531.         vec->ve[i] = mat->me[i][col];
  532.  
  533.     return (vec);
  534. }
  535.  
  536. /* zget_row -- gets a specified row of a matrix and retruns it as a vector */
  537. ZVEC    *zget_row(mat,row,vec)
  538. int    row;
  539. ZMAT    *mat;
  540. ZVEC    *vec;
  541. {
  542.     int    /* i, */ lim;
  543.  
  544.     if ( mat==ZMNULL )
  545.         error(E_NULL,"zget_row");
  546.     if ( row < 0 || row >= mat->m )
  547.         error(E_RANGE,"zget_row");
  548.     if ( vec==ZVNULL || vec->dim<mat->n )
  549.         vec = zv_resize(vec,mat->n);
  550.  
  551.     lim = min(mat->n,vec->dim);
  552.  
  553.     /* for ( i=0; i<mat->n; i++ ) */
  554.     /*     vec->ve[i] = mat->me[row][i]; */
  555.     MEMCOPY(mat->me[row],vec->ve,lim,complex);
  556.  
  557.     return (vec);
  558. }
  559.  
  560. /* zset_col -- sets column of matrix to values given in vec (in situ) */
  561. ZMAT    *zset_col(mat,col,vec)
  562. ZMAT    *mat;
  563. ZVEC    *vec;
  564. int    col;
  565. {
  566.     u_int    i,lim;
  567.  
  568.     if ( mat==ZMNULL || vec==ZVNULL )
  569.         error(E_NULL,"zset_col");
  570.     if ( col < 0 || col >= mat->n )
  571.         error(E_RANGE,"zset_col");
  572.     lim = min(mat->m,vec->dim);
  573.     for ( i=0; i<lim; i++ )
  574.         mat->me[i][col] = vec->ve[i];
  575.  
  576.     return (mat);
  577. }
  578.  
  579. /* zset_row -- sets row of matrix to values given in vec (in situ) */
  580. ZMAT    *zset_row(mat,row,vec)
  581. ZMAT    *mat;
  582. ZVEC    *vec;
  583. int    row;
  584. {
  585.     u_int    /* j, */ lim;
  586.  
  587.     if ( mat==ZMNULL || vec==ZVNULL )
  588.         error(E_NULL,"zset_row");
  589.     if ( row < 0 || row >= mat->m )
  590.         error(E_RANGE,"zset_row");
  591.     lim = min(mat->n,vec->dim);
  592.     /* for ( j=j0; j<lim; j++ ) */
  593.     /*     mat->me[row][j] = vec->ve[j]; */
  594.     MEMCOPY(vec->ve,mat->me[row],lim,complex);
  595.  
  596.     return (mat);
  597. }
  598.  
  599. /* zm_rand -- randomise a complex matrix; uniform in [0,1)+[0,1)*i */
  600. ZMAT    *zm_rand(A)
  601. ZMAT    *A;
  602. {
  603.     int        i;
  604.  
  605.     if ( ! A )
  606.     error(E_NULL,"zm_rand");
  607.  
  608.     for ( i = 0; i < A->m; i++ )
  609.     mrandlist((Real *)(A->me[i]),2*A->n);
  610.  
  611.     return A;
  612. }
  613.