zgivens.c 4.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184
  1. /**************************************************************************
  2. **
  3. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  4. **
  5. ** Meschach Library
  6. **
  7. ** This Meschach Library is provided "as is" without any express
  8. ** or implied warranty of any kind with respect to this software.
  9. ** In particular the authors shall not be liable for any direct,
  10. ** indirect, special, incidental or consequential damages arising
  11. ** in any way from use of the software.
  12. **
  13. ** Everyone is granted permission to copy, modify and redistribute this
  14. ** Meschach Library, provided:
  15. ** 1. All copies contain this copyright notice.
  16. ** 2. All modified copies shall carry a notice stating who
  17. ** made the last modification and the date of such modification.
  18. ** 3. No charge is made for this software or works derived from it.
  19. ** This clause shall not be construed as constraining other software
  20. ** distributed on the same medium as this software, nor is a
  21. ** distribution fee considered a charge.
  22. **
  23. ***************************************************************************/
  24. /*
  25. Givens operations file. Contains routines for calculating and
  26. applying givens rotations for/to vectors and also to matrices by
  27. row and by column.
  28. Complex version.
  29. */
  30. static char rcsid[] = "$Id: ";
  31. #include <stdio.h>
  32. #include <math.h>
  33. #include "zmatrix.h"
  34. #include "zmatrix2.h"
  35. /*
  36. (Complex) Givens rotation matrix:
  37. [ c -s ]
  38. [ s* c ]
  39. Note that c is real and s is complex
  40. */
  41. /* zgivens -- returns c,s parameters for Givens rotation to
  42. eliminate y in the **column** vector [ x y ] */
  43. void zgivens(x,y,c,s)
  44. complex x,y,*s;
  45. Real *c;
  46. {
  47. Real inv_norm, norm;
  48. complex tmp;
  49. /* this is a safe way of computing sqrt(|x|^2+|y|^2) */
  50. tmp.re = zabs(x); tmp.im = zabs(y);
  51. norm = zabs(tmp);
  52. if ( norm == 0.0 )
  53. { *c = 1.0; s->re = s->im = 0.0; } /* identity */
  54. else
  55. {
  56. inv_norm = 1.0 / tmp.re; /* inv_norm = 1/|x| */
  57. x.re *= inv_norm;
  58. x.im *= inv_norm; /* normalise x */
  59. inv_norm = 1.0/norm; /* inv_norm = 1/||[x,y]||2 */
  60. *c = tmp.re * inv_norm;
  61. /* now compute - conj(normalised x).y/||[x,y]||2 */
  62. s->re = - inv_norm*(x.re*y.re + x.im*y.im);
  63. s->im = inv_norm*(x.re*y.im - x.im*y.re);
  64. }
  65. }
  66. /* rot_zvec -- apply Givens rotation to x's i & k components */
  67. ZVEC *rot_zvec(x,i,k,c,s,out)
  68. ZVEC *x,*out;
  69. int i,k;
  70. double c;
  71. complex s;
  72. {
  73. complex temp1, temp2;
  74. if ( x==ZVNULL )
  75. error(E_NULL,"rot_zvec");
  76. if ( i < 0 || i >= x->dim || k < 0 || k >= x->dim )
  77. error(E_RANGE,"rot_zvec");
  78. if ( x != out )
  79. out = zv_copy(x,out);
  80. /* temp1 = c*out->ve[i] - s*out->ve[k]; */
  81. temp1.re = c*out->ve[i].re
  82. - s.re*out->ve[k].re + s.im*out->ve[k].im;
  83. temp1.im = c*out->ve[i].im
  84. - s.re*out->ve[k].im - s.im*out->ve[k].re;
  85. /* temp2 = c*out->ve[k] + zconj(s)*out->ve[i]; */
  86. temp2.re = c*out->ve[k].re
  87. + s.re*out->ve[i].re + s.im*out->ve[i].im;
  88. temp2.im = c*out->ve[k].im
  89. + s.re*out->ve[i].im - s.im*out->ve[i].re;
  90. out->ve[i] = temp1;
  91. out->ve[k] = temp2;
  92. return (out);
  93. }
  94. /* zrot_rows -- premultiply mat by givens rotation described by c,s */
  95. ZMAT *zrot_rows(mat,i,k,c,s,out)
  96. ZMAT *mat,*out;
  97. int i,k;
  98. double c;
  99. complex s;
  100. {
  101. unsigned int j;
  102. complex temp1, temp2;
  103. if ( mat==ZMNULL )
  104. error(E_NULL,"zrot_rows");
  105. if ( i < 0 || i >= mat->m || k < 0 || k >= mat->m )
  106. error(E_RANGE,"zrot_rows");
  107. if ( mat != out )
  108. out = zm_copy(mat,zm_resize(out,mat->m,mat->n));
  109. /* temp1 = c*out->me[i][j] - s*out->me[k][j]; */
  110. for ( j=0; j<mat->n; j++ )
  111. {
  112. /* temp1 = c*out->me[i][j] - s*out->me[k][j]; */
  113. temp1.re = c*out->me[i][j].re
  114. - s.re*out->me[k][j].re + s.im*out->me[k][j].im;
  115. temp1.im = c*out->me[i][j].im
  116. - s.re*out->me[k][j].im - s.im*out->me[k][j].re;
  117. /* temp2 = c*out->me[k][j] + conj(s)*out->me[i][j]; */
  118. temp2.re = c*out->me[k][j].re
  119. + s.re*out->me[i][j].re + s.im*out->me[i][j].im;
  120. temp2.im = c*out->me[k][j].im
  121. + s.re*out->me[i][j].im - s.im*out->me[i][j].re;
  122. out->me[i][j] = temp1;
  123. out->me[k][j] = temp2;
  124. }
  125. return (out);
  126. }
  127. /* zrot_cols -- postmultiply mat by adjoint Givens rotation described by c,s */
  128. ZMAT *zrot_cols(mat,i,k,c,s,out)
  129. ZMAT *mat,*out;
  130. int i,k;
  131. double c;
  132. complex s;
  133. {
  134. unsigned int j;
  135. complex x, y;
  136. if ( mat==ZMNULL )
  137. error(E_NULL,"zrot_cols");
  138. if ( i < 0 || i >= mat->n || k < 0 || k >= mat->n )
  139. error(E_RANGE,"zrot_cols");
  140. if ( mat != out )
  141. out = zm_copy(mat,zm_resize(out,mat->m,mat->n));
  142. for ( j=0; j<mat->m; j++ )
  143. {
  144. x = out->me[j][i]; y = out->me[j][k];
  145. /* out->me[j][i] = c*x - conj(s)*y; */
  146. out->me[j][i].re = c*x.re - s.re*y.re - s.im*y.im;
  147. out->me[j][i].im = c*x.im - s.re*y.im + s.im*y.re;
  148. /* out->me[j][k] = c*y + s*x; */
  149. out->me[j][k].re = c*y.re + s.re*x.re - s.im*x.im;
  150. out->me[j][k].im = c*y.im + s.re*x.im + s.im*x.re;
  151. }
  152. return (out);
  153. }