home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / mesch12a.zip / zhessen.c < prev    next >
C/C++ Source or Header  |  1994-01-13  |  4KB  |  153 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.         File containing routines for determining Hessenberg
  29.     factorisations.
  30.  
  31.     Complex version
  32. */
  33.  
  34. static    char    rcsid[] = "$Id: zhessen.c,v 1.1 1994/01/13 04:27:00 des Exp $";
  35.  
  36. #include    <stdio.h>
  37. #include    "zmatrix.h"
  38. #include        "zmatrix2.h"
  39.  
  40.  
  41. /* zHfactor -- compute Hessenberg factorisation in compact form.
  42.     -- factorisation performed in situ
  43.     -- for details of the compact form see zQRfactor.c and zmatrix2.doc */
  44. ZMAT    *zHfactor(A, diag)
  45. ZMAT    *A;
  46. ZVEC    *diag;
  47. {
  48.     static    ZVEC    *tmp1 = ZVNULL;
  49.     Real    beta;
  50.     int    k, limit;
  51.  
  52.     if ( ! A || ! diag )
  53.         error(E_NULL,"zHfactor");
  54.     if ( diag->dim < A->m - 1 )
  55.         error(E_SIZES,"zHfactor");
  56.     if ( A->m != A->n )
  57.         error(E_SQUARE,"zHfactor");
  58.     limit = A->m - 1;
  59.  
  60.     tmp1 = zv_resize(tmp1,A->m);
  61.     MEM_STAT_REG(tmp1,TYPE_ZVEC);
  62.  
  63.     for ( k = 0; k < limit; k++ )
  64.     {
  65.         zget_col(A,k,tmp1);
  66.         zhhvec(tmp1,k+1,&beta,tmp1,&A->me[k+1][k]);
  67.         diag->ve[k] = tmp1->ve[k+1];
  68.         /* printf("zHfactor: k = %d, beta = %g, tmp1 =\n",k,beta);
  69.         zv_output(tmp1); */
  70.         
  71.         zhhtrcols(A,k+1,k+1,tmp1,beta);
  72.         zhhtrrows(A,0  ,k+1,tmp1,beta);
  73.         /* printf("# at stage k = %d, A =\n",k);    zm_output(A); */
  74.     }
  75.  
  76.     return (A);
  77. }
  78.  
  79. /* zHQunpack -- unpack the compact representation of H and Q of a
  80.     Hessenberg factorisation
  81.     -- if either H or Q is NULL, then it is not unpacked
  82.     -- it can be in situ with HQ == H
  83.     -- returns HQ
  84. */
  85. ZMAT    *zHQunpack(HQ,diag,Q,H)
  86. ZMAT    *HQ, *Q, *H;
  87. ZVEC    *diag;
  88. {
  89.     int    i, j, limit;
  90.     Real    beta, r_ii, tmp_val;
  91.     static    ZVEC    *tmp1 = ZVNULL, *tmp2 = ZVNULL;
  92.  
  93.     if ( HQ==ZMNULL || diag==ZVNULL )
  94.         error(E_NULL,"zHQunpack");
  95.     if ( HQ == Q || H == Q )
  96.         error(E_INSITU,"zHQunpack");
  97.     limit = HQ->m - 1;
  98.     if ( diag->dim < limit )
  99.         error(E_SIZES,"zHQunpack");
  100.     if ( HQ->m != HQ->n )
  101.         error(E_SQUARE,"zHQunpack");
  102.  
  103.  
  104.     if ( Q != ZMNULL )
  105.     {
  106.         Q = zm_resize(Q,HQ->m,HQ->m);
  107.         tmp1 = zv_resize(tmp1,H->m);
  108.         tmp2 = zv_resize(tmp2,H->m);
  109.         MEM_STAT_REG(tmp1,TYPE_ZVEC);
  110.         MEM_STAT_REG(tmp2,TYPE_ZVEC);
  111.         
  112.         for ( i = 0; i < H->m; i++ )
  113.         {
  114.         /* tmp1 = i'th basis vector */
  115.         for ( j = 0; j < H->m; j++ )
  116.             tmp1->ve[j].re = tmp1->ve[j].im = 0.0;
  117.         tmp1->ve[i].re = 1.0;
  118.         
  119.         /* apply H/h transforms in reverse order */
  120.         for ( j = limit-1; j >= 0; j-- )
  121.         {
  122.             zget_col(HQ,j,tmp2);
  123.             r_ii = zabs(tmp2->ve[j+1]);
  124.             tmp2->ve[j+1] = diag->ve[j];
  125.             tmp_val = (r_ii*zabs(diag->ve[j]));
  126.             beta = ( tmp_val == 0.0 ) ? 0.0 : 1.0/tmp_val;
  127.             /* printf("zHQunpack: j = %d, beta = %g, tmp2 =\n",
  128.                j,beta);
  129.             zv_output(tmp2); */
  130.             zhhtrvec(tmp2,beta,j+1,tmp1,tmp1);
  131.         }
  132.         
  133.         /* insert into Q */
  134.         zset_col(Q,i,tmp1);
  135.         }
  136.     }
  137.  
  138.     if ( H != ZMNULL )
  139.     {
  140.         H = zm_copy(HQ,H);
  141.         
  142.         limit = H->m;
  143.         for ( i = 1; i < limit; i++ )
  144.         for ( j = 0; j < i-1; j++ )
  145.             H->me[i][j].re = H->me[i][j].im = 0.0;
  146.     }
  147.  
  148.     return HQ;
  149. }
  150.  
  151.  
  152.  
  153.