home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / mesch12a.zip / zhsehldr.c < prev    next >
C/C++ Source or Header  |  1994-01-13  |  6KB  |  201 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.1 1994/01/13 04:26:43 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.     abs_val = -norm / abs_val;
  71.     newval->re = abs_val*out->ve[i0].re;
  72.     newval->im = abs_val*out->ve[i0].im;
  73.     out->ve[i0].re -= newval->re;
  74.     out->ve[i0].im -= newval->im;
  75.  
  76.     return (out);
  77. }
  78.  
  79. /* zhhtrvec -- apply Householder transformation to vector -- may be in-situ */
  80. ZVEC    *zhhtrvec(hh,beta,i0,in,out)
  81. ZVEC    *hh,*in,*out;    /* hh = Householder vector */
  82. int    i0;
  83. double    beta;
  84. {
  85.     complex    scale, tmp;
  86.     /* u_int    i; */
  87.  
  88.     if ( hh==ZVNULL || in==ZVNULL )
  89.         error(E_NULL,"zhhtrvec");
  90.     if ( in->dim != hh->dim )
  91.         error(E_SIZES,"zhhtrvec");
  92.     if ( i0 < 0 || i0 > in->dim )
  93.         error(E_BOUNDS,"zhhvec");
  94.  
  95.     tmp = _zin_prod(hh,in,i0,Z_CONJ);
  96.     scale.re = -beta*tmp.re;
  97.     scale.im = -beta*tmp.im;
  98.     out = zv_copy(in,out);
  99.     __zmltadd__(&(out->ve[i0]),&(hh->ve[i0]),scale,
  100.             (int)(in->dim-i0),Z_NOCONJ);
  101.     /************************************************************
  102.     for ( i=i0; i<in->dim; i++ )
  103.         out->ve[i] = in->ve[i] - scale*hh->ve[i];
  104.     ************************************************************/
  105.  
  106.     return (out);
  107. }
  108.  
  109. /* zhhtrrows -- transform a matrix by a Householder vector by rows
  110.     starting at row i0 from column j0 -- in-situ */
  111. ZMAT    *zhhtrrows(M,i0,j0,hh,beta)
  112. ZMAT    *M;
  113. int    i0, j0;
  114. ZVEC    *hh;
  115. double    beta;
  116. {
  117.     complex    ip, scale;
  118.     int    i /*, j */;
  119.  
  120.     if ( M==ZMNULL || hh==ZVNULL )
  121.         error(E_NULL,"zhhtrrows");
  122.     if ( M->n != hh->dim )
  123.         error(E_RANGE,"zhhtrrows");
  124.     if ( i0 < 0 || i0 > M->m || j0 < 0 || j0 > M->n )
  125.         error(E_BOUNDS,"zhhtrrows");
  126.  
  127.     if ( beta == 0.0 )    return (M);
  128.  
  129.     /* for each row ... */
  130.     for ( i = i0; i < M->m; i++ )
  131.     {    /* compute inner product */
  132.         ip = __zip__(&(M->me[i][j0]),&(hh->ve[j0]),
  133.                  (int)(M->n-j0),Z_NOCONJ);
  134.         /**************************************************
  135.         ip = 0.0;
  136.         for ( j = j0; j < M->n; j++ )
  137.             ip += M->me[i][j]*hh->ve[j];
  138.         **************************************************/
  139.         scale.re = -beta*ip.re;
  140.         scale.im = -beta*ip.im;
  141.         /* if ( scale == 0.0 ) */
  142.         if ( is_zero(scale) )
  143.             continue;
  144.  
  145.         /* do operation */
  146.         __zmltadd__(&(M->me[i][j0]),&(hh->ve[j0]),scale,
  147.                 (int)(M->n-j0),Z_CONJ);
  148.         /**************************************************
  149.         for ( j = j0; j < M->n; j++ )
  150.             M->me[i][j] -= scale*hh->ve[j];
  151.         **************************************************/
  152.     }
  153.  
  154.     return (M);
  155. }
  156.  
  157.  
  158. /* zhhtrcols -- transform a matrix by a Householder vector by columns
  159.     starting at row i0 from column j0 -- in-situ */
  160. ZMAT    *zhhtrcols(M,i0,j0,hh,beta)
  161. ZMAT    *M;
  162. int    i0, j0;
  163. ZVEC    *hh;
  164. double    beta;
  165. {
  166.     /* Real    ip, scale; */
  167.     complex    scale;
  168.     int    i /*, k */;
  169.     static    ZVEC    *w = ZVNULL;
  170.  
  171.     if ( M==ZMNULL || hh==ZVNULL )
  172.         error(E_NULL,"zhhtrcols");
  173.     if ( M->m != hh->dim )
  174.         error(E_SIZES,"zhhtrcols");
  175.     if ( i0 < 0 || i0 > M->m || j0 < 0 || j0 > M->n )
  176.         error(E_BOUNDS,"zhhtrcols");
  177.  
  178.     if ( beta == 0.0 )    return (M);
  179.  
  180.     w = zv_resize(w,M->n);
  181.     MEM_STAT_REG(w,TYPE_ZVEC);
  182.     zv_zero(w);
  183.  
  184.     for ( i = i0; i < M->m; i++ )
  185.         /* if ( hh->ve[i] != 0.0 ) */
  186.         if ( ! is_zero(hh->ve[i]) )
  187.         __zmltadd__(&(w->ve[j0]),&(M->me[i][j0]),hh->ve[i],
  188.                 (int)(M->n-j0),Z_CONJ);
  189.     for ( i = i0; i < M->m; i++ )
  190.         /* if ( hh->ve[i] != 0.0 ) */
  191.         if ( ! is_zero(hh->ve[i]) )
  192.         {
  193.         scale.re = -beta*hh->ve[i].re;
  194.         scale.im = -beta*hh->ve[i].im;
  195.         __zmltadd__(&(M->me[i][j0]),&(w->ve[j0]),scale,
  196.                 (int)(M->n-j0),Z_CONJ);
  197.         }
  198.     return (M);
  199. }
  200.  
  201.