home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / mesch12a.zip / hessen.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. /*
  29.         File containing routines for determining Hessenberg
  30.     factorisations.
  31. */
  32.  
  33. static    char    rcsid[] = "$Id: hessen.c,v 1.2 1994/01/13 05:36:24 des Exp $";
  34.  
  35. #include    <stdio.h>
  36. #include    "matrix.h"
  37. #include        "matrix2.h"
  38.  
  39.  
  40.  
  41. /* Hfactor -- compute Hessenberg factorisation in compact form.
  42.     -- factorisation performed in situ
  43.     -- for details of the compact form see QRfactor.c and matrix2.doc */
  44. MAT    *Hfactor(A, diag, beta)
  45. MAT    *A;
  46. VEC    *diag, *beta;
  47. {
  48.     static    VEC    *tmp1 = VNULL;
  49.     int    k, limit;
  50.  
  51.     if ( ! A || ! diag || ! beta )
  52.         error(E_NULL,"Hfactor");
  53.     if ( diag->dim < A->m - 1 || beta->dim < A->m - 1 )
  54.         error(E_SIZES,"Hfactor");
  55.     if ( A->m != A->n )
  56.         error(E_SQUARE,"Hfactor");
  57.     limit = A->m - 1;
  58.  
  59.     tmp1 = v_resize(tmp1,A->m);
  60.     MEM_STAT_REG(tmp1,TYPE_VEC);
  61.  
  62.     for ( k = 0; k < limit; k++ )
  63.     {
  64.         get_col(A,(u_int)k,tmp1);
  65.         /* printf("the %d'th column = ");    v_output(tmp1); */
  66.         hhvec(tmp1,k+1,&beta->ve[k],tmp1,&A->me[k+1][k]);
  67.         /* diag->ve[k] = tmp1->ve[k+1]; */
  68.         v_set_val(diag,k,v_entry(tmp1,k+1));
  69.         /* printf("H/h vector = ");    v_output(tmp1); */
  70.         /* printf("from the %d'th entry\n",k+1); */
  71.         /* printf("beta = %g\n",beta->ve[k]); */
  72.  
  73.         /* hhtrcols(A,k+1,k+1,tmp1,beta->ve[k]); */
  74.         /* hhtrrows(A,0  ,k+1,tmp1,beta->ve[k]); */
  75.         hhtrcols(A,k+1,k+1,tmp1,v_entry(beta,k));
  76.         hhtrrows(A,0  ,k+1,tmp1,v_entry(beta,k));
  77.         /* printf("A = ");        m_output(A); */
  78.     }
  79.  
  80.     return (A);
  81. }
  82.  
  83. /* makeHQ -- construct the Hessenberg orthogonalising matrix Q;
  84.     -- i.e. Hess M = Q.M.Q'    */
  85. MAT    *makeHQ(H, diag, beta, Qout)
  86. MAT    *H, *Qout;
  87. VEC    *diag, *beta;
  88. {
  89.     int    i, j, limit;
  90.     static    VEC    *tmp1 = VNULL, *tmp2 = VNULL;
  91.  
  92.     if ( H==(MAT *)NULL || diag==(VEC *)NULL || beta==(VEC *)NULL )
  93.         error(E_NULL,"makeHQ");
  94.     limit = H->m - 1;
  95.     if ( diag->dim < limit || beta->dim < limit )
  96.         error(E_SIZES,"makeHQ");
  97.     if ( H->m != H->n )
  98.         error(E_SQUARE,"makeHQ");
  99.     Qout = m_resize(Qout,H->m,H->m);
  100.  
  101.     tmp1 = v_resize(tmp1,H->m);
  102.     tmp2 = v_resize(tmp2,H->m);
  103.     MEM_STAT_REG(tmp1,TYPE_VEC);
  104.     MEM_STAT_REG(tmp2,TYPE_VEC);
  105.  
  106.     for ( i = 0; i < H->m; i++ )
  107.     {
  108.         /* tmp1 = i'th basis vector */
  109.         for ( j = 0; j < H->m; j++ )
  110.             /* tmp1->ve[j] = 0.0; */
  111.             v_set_val(tmp1,j,0.0);
  112.         /* tmp1->ve[i] = 1.0; */
  113.         v_set_val(tmp1,i,1.0);
  114.  
  115.         /* apply H/h transforms in reverse order */
  116.         for ( j = limit-1; j >= 0; j-- )
  117.         {
  118.             get_col(H,(u_int)j,tmp2);
  119.             /* tmp2->ve[j+1] = diag->ve[j]; */
  120.             v_set_val(tmp2,j+1,v_entry(diag,j));
  121.             hhtrvec(tmp2,beta->ve[j],j+1,tmp1,tmp1);
  122.         }
  123.  
  124.         /* insert into Qout */
  125.         set_col(Qout,(u_int)i,tmp1);
  126.     }
  127.  
  128.     return (Qout);
  129. }
  130.  
  131. /* makeH -- construct actual Hessenberg matrix */
  132. MAT    *makeH(H,Hout)
  133. MAT    *H, *Hout;
  134. {
  135.     int    i, j, limit;
  136.  
  137.     if ( H==(MAT *)NULL )
  138.         error(E_NULL,"makeH");
  139.     if ( H->m != H->n )
  140.         error(E_SQUARE,"makeH");
  141.     Hout = m_resize(Hout,H->m,H->m);
  142.     Hout = m_copy(H,Hout);
  143.  
  144.     limit = H->m;
  145.     for ( i = 1; i < limit; i++ )
  146.         for ( j = 0; j < i-1; j++ )
  147.             /* Hout->me[i][j] = 0.0;*/
  148.             m_set_val(Hout,i,j,0.0);
  149.  
  150.     return (Hout);
  151. }
  152.  
  153.