home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / mesch12a.zip / givens.c < prev    next >
C/C++ Source or Header  |  1994-01-13  |  4KB  |  139 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. /*
  29.         Files for matrix computations
  30.  
  31.     Givens operations file. Contains routines for calculating and
  32.     applying givens rotations for/to vectors and also to matrices by
  33.     row and by column.
  34. */
  35.  
  36. /* givens.c 1.2 11/25/87 */
  37. static    char    rcsid[] = "$Id: givens.c,v 1.2 1994/01/13 05:39:42 des Exp $";
  38.  
  39. #include    <stdio.h>
  40. #include    <math.h>
  41. #include    "matrix.h"
  42. #include        "matrix2.h"
  43.  
  44. /* givens -- returns c,s parameters for Givens rotation to
  45.         eliminate y in the vector [ x y ]' */
  46. void    givens(x,y,c,s)
  47. double  x,y;
  48. Real    *c,*s;
  49. {
  50.     Real    norm;
  51.  
  52.     norm = sqrt(x*x+y*y);
  53.     if ( norm == 0.0 )
  54.     {    *c = 1.0;    *s = 0.0;    }    /* identity */
  55.     else
  56.     {    *c = x/norm;    *s = y/norm;    }
  57. }
  58.  
  59. /* rot_vec -- apply Givens rotation to x's i & k components */
  60. VEC    *rot_vec(x,i,k,c,s,out)
  61. VEC    *x,*out;
  62. u_int    i,k;
  63. double    c,s;
  64. {
  65.     Real    temp;
  66.  
  67.     if ( x==VNULL )
  68.         error(E_NULL,"rot_vec");
  69.     if ( i >= x->dim || k >= x->dim )
  70.         error(E_RANGE,"rot_vec");
  71.     out = v_copy(x,out);
  72.  
  73.     /* temp = c*out->ve[i] + s*out->ve[k]; */
  74.     temp = c*v_entry(out,i) + s*v_entry(out,k);
  75.     /* out->ve[k] = -s*out->ve[i] + c*out->ve[k]; */
  76.     v_set_val(out,k,-s*v_entry(out,i)+c*v_entry(out,k));
  77.     /* out->ve[i] = temp; */
  78.     v_set_val(out,i,temp);
  79.  
  80.     return (out);
  81. }
  82.  
  83. /* rot_rows -- premultiply mat by givens rotation described by c,s */
  84. MAT    *rot_rows(mat,i,k,c,s,out)
  85. MAT    *mat,*out;
  86. u_int    i,k;
  87. double    c,s;
  88. {
  89.     u_int    j;
  90.     Real    temp;
  91.  
  92.     if ( mat==(MAT *)NULL )
  93.         error(E_NULL,"rot_rows");
  94.     if ( i >= mat->m || k >= mat->m )
  95.         error(E_RANGE,"rot_rows");
  96.     out = m_copy(mat,out);
  97.  
  98.     for ( j=0; j<mat->n; j++ )
  99.     {
  100.         /* temp = c*out->me[i][j] + s*out->me[k][j]; */
  101.         temp = c*m_entry(out,i,j) + s*m_entry(out,k,j);
  102.         /* out->me[k][j] = -s*out->me[i][j] + c*out->me[k][j]; */
  103.         m_set_val(out,k,j, -s*m_entry(out,i,j) + c*m_entry(out,k,j));
  104.         /* out->me[i][j] = temp; */
  105.         m_set_val(out,i,j, temp);
  106.     }
  107.  
  108.     return (out);
  109. }
  110.  
  111. /* rot_cols -- postmultiply mat by givens rotation described by c,s */
  112. MAT    *rot_cols(mat,i,k,c,s,out)
  113. MAT    *mat,*out;
  114. u_int    i,k;
  115. double    c,s;
  116. {
  117.     u_int    j;
  118.     Real    temp;
  119.  
  120.     if ( mat==(MAT *)NULL )
  121.         error(E_NULL,"rot_cols");
  122.     if ( i >= mat->n || k >= mat->n )
  123.         error(E_RANGE,"rot_cols");
  124.     out = m_copy(mat,out);
  125.  
  126.     for ( j=0; j<mat->m; j++ )
  127.     {
  128.         /* temp = c*out->me[j][i] + s*out->me[j][k]; */
  129.         temp = c*m_entry(out,j,i) + s*m_entry(out,j,k);
  130.         /* out->me[j][k] = -s*out->me[j][i] + c*out->me[j][k]; */
  131.         m_set_val(out,j,k, -s*m_entry(out,j,i) + c*m_entry(out,j,k));
  132.         /* out->me[j][i] = temp; */
  133.         m_set_val(out,j,i,temp);
  134.     }
  135.  
  136.     return (out);
  137. }
  138.  
  139.