home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / mesch12a.zip / zgivens.c < prev    next >
C/C++ Source or Header  |  1994-01-13  |  5KB  |  182 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.     Givens operations file. Contains routines for calculating and
  29.     applying givens rotations for/to vectors and also to matrices by
  30.     row and by column.
  31.  
  32.     Complex version.
  33. */
  34.  
  35. static    char    rcsid[] = "$Id: ";
  36.  
  37. #include    <stdio.h>
  38. #include    <math.h>
  39. #include    "zmatrix.h"
  40. #include        "zmatrix2.h"
  41.  
  42. /*
  43.     (Complex) Givens rotation matrix:
  44.         [ c   -s ]
  45.         [ s*   c ]
  46.     Note that c is real and s is complex
  47. */
  48.  
  49. /* zgivens -- returns c,s parameters for Givens rotation to
  50.         eliminate y in the **column** vector [ x y ] */
  51. void    zgivens(x,y,c,s)
  52. complex    x,y,*s;
  53. Real    *c;
  54. {
  55.     Real    inv_norm, norm;
  56.     complex    old_x, tmp;
  57.  
  58.     /* this is a safe way of computing sqrt(|x|^2+|y|^2) */
  59.     tmp.re = zabs(x);    tmp.im = zabs(y);
  60.     norm = zabs(tmp);
  61.  
  62.     if ( norm == 0.0 )
  63.     {    *c = 1.0;    s->re = s->im = 0.0;    } /* identity */
  64.     else
  65.     {
  66.         inv_norm = 1.0 / tmp.re;    /* inv_norm = 1/|x| */
  67.         old_x = x;
  68.         x.re *= inv_norm;
  69.         x.im *= inv_norm;        /* normalise x */
  70.         inv_norm = 1.0/norm;        /* inv_norm = 1/||[x,y]||2 */
  71.         *c = tmp.re * inv_norm;
  72.         /* now compute - conj(normalised x).y/||[x,y]||2 */
  73.         s->re = - inv_norm*(x.re*y.re + x.im*y.im);
  74.         s->im =   inv_norm*(x.re*y.im - x.im*y.re);
  75.     }
  76. }
  77.  
  78. /* rot_zvec -- apply Givens rotation to x's i & k components */
  79. ZVEC    *rot_zvec(x,i,k,c,s,out)
  80. ZVEC    *x,*out;
  81. int    i,k;
  82. double    c;
  83. complex    s;
  84. {
  85.  
  86.     complex    temp1, temp2;
  87.  
  88.     if ( x==ZVNULL )
  89.         error(E_NULL,"rot_zvec");
  90.     if ( i < 0 || i >= x->dim || k < 0 || k >= x->dim )
  91.         error(E_RANGE,"rot_zvec");
  92.     if ( x != out )
  93.         out = zv_copy(x,out);
  94.  
  95.     /* temp1 = c*out->ve[i] - s*out->ve[k]; */
  96.     temp1.re = c*out->ve[i].re
  97.         - s.re*out->ve[k].re + s.im*out->ve[k].im;
  98.     temp1.im = c*out->ve[i].im
  99.         - s.re*out->ve[k].im - s.im*out->ve[k].re;
  100.  
  101.     /* temp2 = c*out->ve[k] + zconj(s)*out->ve[i]; */
  102.     temp2.re = c*out->ve[k].re
  103.         + s.re*out->ve[i].re + s.im*out->ve[i].im;
  104.     temp2.im = c*out->ve[k].im
  105.         + s.re*out->ve[i].im - s.im*out->ve[i].re;
  106.  
  107.     out->ve[i] = temp1;
  108.     out->ve[k] = temp2;
  109.  
  110.     return (out);
  111. }
  112.  
  113. /* zrot_rows -- premultiply mat by givens rotation described by c,s */
  114. ZMAT    *zrot_rows(mat,i,k,c,s,out)
  115. ZMAT    *mat,*out;
  116. int    i,k;
  117. double    c;
  118. complex    s;
  119. {
  120.     u_int    j;
  121.     complex    temp1, temp2;
  122.  
  123.     if ( mat==ZMNULL )
  124.         error(E_NULL,"zrot_rows");
  125.     if ( i < 0 || i >= mat->m || k < 0 || k >= mat->m )
  126.         error(E_RANGE,"zrot_rows");
  127.     out = zm_copy(mat,out);
  128.  
  129.     /* temp1 = c*out->me[i][j] - s*out->me[k][j]; */
  130.     for ( j=0; j<mat->n; j++ )
  131.     {
  132.         /* temp1 = c*out->me[i][j] - s*out->me[k][j]; */
  133.         temp1.re = c*out->me[i][j].re
  134.         - s.re*out->me[k][j].re + s.im*out->me[k][j].im;
  135.         temp1.im = c*out->me[i][j].im
  136.         - s.re*out->me[k][j].im - s.im*out->me[k][j].re;
  137.         
  138.         /* temp2 = c*out->me[k][j] + conj(s)*out->me[i][j]; */
  139.         temp2.re = c*out->me[k][j].re
  140.         + s.re*out->me[i][j].re + s.im*out->me[i][j].im;
  141.         temp2.im = c*out->me[k][j].im
  142.         + s.re*out->me[i][j].im - s.im*out->me[i][j].re;
  143.         
  144.         out->me[i][j] = temp1;
  145.         out->me[k][j] = temp2;
  146.     }
  147.  
  148.     return (out);
  149. }
  150.  
  151. /* zrot_cols -- postmultiply mat by adjoint Givens rotation described by c,s */
  152. ZMAT    *zrot_cols(mat,i,k,c,s,out)
  153. ZMAT    *mat,*out;
  154. int    i,k;
  155. double    c;
  156. complex    s;
  157. {
  158.     u_int    j;
  159.     complex    x, y;
  160.  
  161.     if ( mat==ZMNULL )
  162.         error(E_NULL,"zrot_cols");
  163.     if ( i < 0 || i >= mat->n || k < 0 || k >= mat->n )
  164.         error(E_RANGE,"zrot_cols");
  165.     out = zm_copy(mat,out);
  166.  
  167.     for ( j=0; j<mat->m; j++ )
  168.     {
  169.         x = out->me[j][i];    y = out->me[j][k];
  170.         /* out->me[j][i] = c*x - conj(s)*y; */
  171.         out->me[j][i].re = c*x.re - s.re*y.re - s.im*y.im;
  172.         out->me[j][i].im = c*x.im - s.re*y.im + s.im*y.re;
  173.         
  174.         /* out->me[j][k] = c*y + s*x; */
  175.         out->me[j][k].re = c*y.re + s.re*x.re - s.im*x.im;
  176.         out->me[j][k].im = c*y.im + s.re*x.im + s.im*x.re;
  177.     }
  178.  
  179.     return (out);
  180. }
  181.  
  182.