home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 1 / RISC_DISC_1.iso / pd_share / code / meschach / !Meschach / c / zhsehldr < prev    next >
Encoding:
Text File  |  1994-04-06  |  5.4 KB  |  209 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.         Files for matrix computations
  29.  
  30.     Householder transformation file. Contains routines for calculating
  31.     householder transformations, applying them to vectors and matrices
  32.     by both row & column.
  33.  
  34.     Complex version
  35. */
  36.  
  37. static    char    rcsid[] = "$Id: zhsehldr.c,v 1.2 1994/04/07 01:43:47 des Exp $";
  38.  
  39. #include    <stdio.h>
  40. #include    <math.h>
  41. #include    "zmatrix.h"
  42. #include        "zmatrix2.h"
  43.  
  44. #define    is_zero(z)    ((z).re == 0.0 && (z).im == 0.0)
  45.  
  46. /* zhhvec -- calulates Householder vector to eliminate all entries after the
  47.     i0 entry of the vector vec. It is returned as out. May be in-situ */
  48. ZVEC    *zhhvec(vec,i0,beta,out,newval)
  49. ZVEC    *vec,*out;
  50. int    i0;
  51. Real    *beta;
  52. complex    *newval;
  53. {
  54.     complex    tmp;
  55.     Real    norm, abs_val;
  56.  
  57.     if ( i0 < 0 || i0 >= vec->dim )
  58.         error(E_BOUNDS,"zhhvec");
  59.     out = _zv_copy(vec,out,i0);
  60.     tmp = _zin_prod(out,out,i0,Z_CONJ);
  61.     if ( tmp.re <= 0.0 )
  62.     {
  63.         *beta = 0.0;
  64.         *newval = out->ve[i0];
  65.         return (out);
  66.     }
  67.     norm = sqrt(tmp.re);
  68.     abs_val = zabs(out->ve[i0]);
  69.     *beta = 1.0/(norm * (norm+abs_val));
  70.     if ( abs_val == 0.0 )
  71.     {
  72.       newval->re = norm;
  73.       newval->im = 0.0;
  74.     }
  75.     else
  76.     { 
  77.       abs_val = -norm / abs_val;
  78.       newval->re = abs_val*out->ve[i0].re;
  79.       newval->im = abs_val*out->ve[i0].im;
  80.     }    abs_val = -norm / abs_val;
  81.     out->ve[i0].re -= newval->re;
  82.     out->ve[i0].im -= newval->im;
  83.  
  84.     return (out);
  85. }
  86.  
  87. /* zhhtrvec -- apply Householder transformation to vector -- may be in-situ */
  88. ZVEC    *zhhtrvec(hh,beta,i0,in,out)
  89. ZVEC    *hh,*in,*out;    /* hh = Householder vector */
  90. int    i0;
  91. double    beta;
  92. {
  93.     complex    scale, tmp;
  94.     /* u_int    i; */
  95.  
  96.     if ( hh==ZVNULL || in==ZVNULL )
  97.         error(E_NULL,"zhhtrvec");
  98.     if ( in->dim != hh->dim )
  99.         error(E_SIZES,"zhhtrvec");
  100.     if ( i0 < 0 || i0 > in->dim )
  101.         error(E_BOUNDS,"zhhvec");
  102.  
  103.     tmp = _zin_prod(hh,in,i0,Z_CONJ);
  104.     scale.re = -beta*tmp.re;
  105.     scale.im = -beta*tmp.im;
  106.     out = zv_copy(in,out);
  107.     __zmltadd__(&(out->ve[i0]),&(hh->ve[i0]),scale,
  108.             (int)(in->dim-i0),Z_NOCONJ);
  109.     /************************************************************
  110.     for ( i=i0; i<in->dim; i++ )
  111.         out->ve[i] = in->ve[i] - scale*hh->ve[i];
  112.     ************************************************************/
  113.  
  114.     return (out);
  115. }
  116.  
  117. /* zhhtrrows -- transform a matrix by a Householder vector by rows
  118.     starting at row i0 from column j0 -- in-situ */
  119. ZMAT    *zhhtrrows(M,i0,j0,hh,beta)
  120. ZMAT    *M;
  121. int    i0, j0;
  122. ZVEC    *hh;
  123. double    beta;
  124. {
  125.     complex    ip, scale;
  126.     int    i /*, j */;
  127.  
  128.     if ( M==ZMNULL || hh==ZVNULL )
  129.         error(E_NULL,"zhhtrrows");
  130.     if ( M->n != hh->dim )
  131.         error(E_RANGE,"zhhtrrows");
  132.     if ( i0 < 0 || i0 > M->m || j0 < 0 || j0 > M->n )
  133.         error(E_BOUNDS,"zhhtrrows");
  134.  
  135.     if ( beta == 0.0 )    return (M);
  136.  
  137.     /* for each row ... */
  138.     for ( i = i0; i < M->m; i++ )
  139.     {    /* compute inner product */
  140.         ip = __zip__(&(M->me[i][j0]),&(hh->ve[j0]),
  141.                  (int)(M->n-j0),Z_NOCONJ);
  142.         /**************************************************
  143.         ip = 0.0;
  144.         for ( j = j0; j < M->n; j++ )
  145.             ip += M->me[i][j]*hh->ve[j];
  146.         **************************************************/
  147.         scale.re = -beta*ip.re;
  148.         scale.im = -beta*ip.im;
  149.         /* if ( scale == 0.0 ) */
  150.         if ( is_zero(scale) )
  151.             continue;
  152.  
  153.         /* do operation */
  154.         __zmltadd__(&(M->me[i][j0]),&(hh->ve[j0]),scale,
  155.                 (int)(M->n-j0),Z_CONJ);
  156.         /**************************************************
  157.         for ( j = j0; j < M->n; j++ )
  158.             M->me[i][j] -= scale*hh->ve[j];
  159.         **************************************************/
  160.     }
  161.  
  162.     return (M);
  163. }
  164.  
  165.  
  166. /* zhhtrcols -- transform a matrix by a Householder vector by columns
  167.     starting at row i0 from column j0 -- in-situ */
  168. ZMAT    *zhhtrcols(M,i0,j0,hh,beta)
  169. ZMAT    *M;
  170. int    i0, j0;
  171. ZVEC    *hh;
  172. double    beta;
  173. {
  174.     /* Real    ip, scale; */
  175.     complex    scale;
  176.     int    i /*, k */;
  177.     static    ZVEC    *w = ZVNULL;
  178.  
  179.     if ( M==ZMNULL || hh==ZVNULL )
  180.         error(E_NULL,"zhhtrcols");
  181.     if ( M->m != hh->dim )
  182.         error(E_SIZES,"zhhtrcols");
  183.     if ( i0 < 0 || i0 > M->m || j0 < 0 || j0 > M->n )
  184.         error(E_BOUNDS,"zhhtrcols");
  185.  
  186.     if ( beta == 0.0 )    return (M);
  187.  
  188.     w = zv_resize(w,M->n);
  189.     MEM_STAT_REG(w,TYPE_ZVEC);
  190.     zv_zero(w);
  191.  
  192.     for ( i = i0; i < M->m; i++ )
  193.         /* if ( hh->ve[i] != 0.0 ) */
  194.         if ( ! is_zero(hh->ve[i]) )
  195.         __zmltadd__(&(w->ve[j0]),&(M->me[i][j0]),hh->ve[i],
  196.                 (int)(M->n-j0),Z_CONJ);
  197.     for ( i = i0; i < M->m; i++ )
  198.         /* if ( hh->ve[i] != 0.0 ) */
  199.         if ( ! is_zero(hh->ve[i]) )
  200.         {
  201.         scale.re = -beta*hh->ve[i].re;
  202.         scale.im = -beta*hh->ve[i].im;
  203.         __zmltadd__(&(M->me[i][j0]),&(w->ve[j0]),scale,
  204.                 (int)(M->n-j0),Z_CONJ);
  205.         }
  206.     return (M);
  207. }
  208.  
  209.