home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 1 / ARM_CLUB_CD.iso / contents / apps / clib / progs / meschach / !Meschach / c / machine < prev    next >
Encoding:
Text File  |  1994-08-03  |  4.0 KB  |  163 lines
  1. tains routines for calculating
  2.     householder transformations, applying them to vectors and matrices
  3.     by both row & column.
  4. */
  5.  
  6. /* hsehldr.c 1.3 10/8/87 */
  7. static    char    rcsid[] = "$Id: hsehldr.c,v 1.2 1994/01/13 05:36:29 des Exp $";
  8.  
  9. #include    <stdio.h>
  10. #include    <math.h>
  11. #include    "matrix.h"
  12. #include        "matrix2.h"
  13.  
  14.  
  15. /* hhvec -- calulates Householder vector to eliminate all entries after the
  16.     i0 entry of the vector vec. It is returned as out. May be in-situ */
  17. VEC    *hhvec(vec,i0,beta,out,newval)
  18. VEC    *vec,*out;
  19. u_int    i0;
  20. Real    *beta,*newval;
  21. {
  22.     Real    norm;
  23.  
  24.     out = _v_copy(vec,out,i0);
  25.     norm = sqrt(_in_prod(out,out,i0));
  26.     if ( norm <= 0.0 )
  27.     {
  28.         *beta = 0.0;
  29.         return (out);
  30.     }
  31.     *beta = 1.0/(norm * (norm+fabs(out->ve[i0])));
  32.     if ( out->ve[i0] > 0.0 )
  33.         *newval = -norm;
  34.     else
  35.         *newval = norm;
  36.     out->ve[i0] -= *newval;
  37.  
  38.     return (out);
  39. }
  40.  
  41. /* hhtrvec -- apply Householder transformation to vector -- may be in-situ */
  42. VEC    *hhtrvec(hh,beta,i0,in,out)
  43. VEC    *hh,*in,*out;    /* hh = Householder vector */
  44. u_int    i0;
  45. double    beta;
  46. {
  47.     Real    scale;
  48.     /* u_int    i; */
  49.  
  50.     if ( hh==(VEC *)NULL || in==(VEC *)NULL )
  51.         error(E_NULL,"hhtrvec");
  52.     if ( in->dim != hh->dim )
  53.         error(E_SIZES,"hhtrvec");
  54.     if ( i0 > in->dim )
  55.         error(E_BOUNDS,"hhtrvec");
  56.  
  57.     scale = beta*_in_prod(hh,in,i0);
  58.     out = v_copy(in,out);
  59.     __mltadd__(&(out->ve[i0]),&(hh->ve[i0]),-scale,(int)(in->dim-i0));
  60.     /************************************************************
  61.     for ( i=i0; i<in->dim; i++ )
  62.         out->ve[i] = in->ve[i] - scale*hh->ve[i];
  63.     ************************************************************/
  64.  
  65.     return (out);
  66. }
  67.  
  68. /* hhtrrows -- transform a matrix by a Householder vector by rows
  69.     starting at row i0 from column j0 -- in-situ */
  70. MAT    *hhtrrows(M,i0,j0,hh,beta)
  71. MAT    *M;
  72. u_int    i0, j0;
  73. VEC    *hh;
  74. double    beta;
  75. {
  76.     Real    ip, scale;
  77.     int    i /*, j */;
  78.  
  79.     if ( M==(MAT *)NULL || hh==(VEC *)NULL )
  80.         error(E_NULL,"hhtrrows");
  81.     if ( M->n != hh->dim )
  82.         error(E_RANGE,"hhtrrows");
  83.     if ( i0 > M->m || j0 > M->n )
  84.         error(E_BOUNDS,"hhtrrows");
  85.  
  86.     if ( beta == 0.0 )    return (M);
  87.  
  88.     /* for each row ... */
  89.     for ( i = i0; i < M->m; i++ )
  90.     {    /* compute inner product */
  91.         ip = __ip__(&(M->me[i][j0]),&(hh->ve[j0]),(int)(M->n-j0));
  92.         /**************************************************
  93.         ip = 0.0;
  94.         for ( j = j0; j < M->n; j++ )
  95.             ip += M->me[i][j]*hh->ve[j];
  96.         **************************************************/
  97.         scale = beta*ip;
  98.         if ( scale == 0.0 )
  99.             continue;
  100.  
  101.         /* do operation */
  102.         __mltadd__(&(M->me[i][j0]),&(hh->ve[j0]),-scale,
  103.                             (int)(M->n-j0));
  104.         /**************************************************
  105.         for ( j = j0; j < M->n; j++ )
  106.             M->me[i][j] -= scale*hh->ve[j];
  107.         **************************************************/
  108.     }
  109.  
  110.     return (M);
  111. }
  112.  
  113.  
  114. /* hhtrcols -- transform a matrix by a Householder vector by columns
  115.     starting at row i0 from column j0 -- in-situ */
  116. MAT    *hhtrcols(M,i0,j0,hh,beta)
  117. MAT    *M;
  118. u_int    i0, j0;
  119. VEC    *hh;
  120. double    beta;
  121. {
  122.     /* Real    ip, scale; */
  123.     int    i /*, k */;
  124.     static    VEC    *w = VNULL;
  125.  
  126.     if ( M==(MAT *)NULL || hh==(VEC *)NULL )
  127.         error(E_NULL,"hhtrcols");
  128.     if ( M->m != hh->dim )
  129.         error(E_SIZES,"hhtrcols");
  130.     if ( i0 > M->m || j0 > M->n )
  131.         error(E_BOUNDS,"hhtrcols");
  132.  
  133.     if ( beta == 0.0 )    return (M);
  134.  
  135.     w = v_resize(w,M->n);
  136.     MEM_STAT_REG(w,TYPE_VEC);
  137.     v_zero(w);
  138.  
  139.     for ( i = i0; i < M->m; i++ )
  140.         if ( hh->ve[i] != 0.0 )
  141.         __mltadd__(&(w->ve[j0]),&(M->me[i][j0]),hh->ve[i],
  142.                             (int)(M->n-j0));
  143.     for ( i = i0; i < M->m; i++ )
  144.         if ( hh->ve[i] != 0.0 )
  145.         __mltadd__(&(M->me[i][j0]),&(w->ve[j0]),-beta*hh->ve[i],
  146.                             (int)(M->n-j0));
  147.     return (M);
  148. }
  149.  
  150. FileDataŵinitEÿÿÿè%@Ðë
  151. /**************************************************************************
  152. **
  153. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  154. **
  155. **                 Meschach Library
  156. ** 
  157. ** This Meschach Library is provided "as is" without any express 
  158. ** or implied warranty of any kind with respect to this software. 
  159. ** In particular the authors shall not be liable for any direct, 
  160. ** indirect, special, incidental or consequential damages arising 
  161. ** in any way from use of the software.
  162. ** 
  163. ** Everyone is granted permission to copy, modify and redistribut