home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 1 / ARM_CLUB_CD.iso / contents / apps / clib / progs / meschach / !Meschach / c / tutorial < prev    next >
Encoding:
Text File  |  1994-05-19  |  7.6 KB  |  291 lines
  1. & ( U->m != U->n ) ) || ( V && ( V->m != V->n ) ) )
  2.         error(E_SQUARE,"bifactor");
  3.     if ( ( U && U->m != A->m ) || ( V && V->m != A->n ) )
  4.         error(E_SIZES,"bifactor");
  5.     tmp1 = v_resize(tmp1,A->m);
  6.     tmp2 = v_resize(tmp2,A->n);
  7.     MEM_STAT_REG(tmp1,TYPE_VEC);
  8.     MEM_STAT_REG(tmp2,TYPE_VEC);
  9.  
  10.     if ( A->m >= A->n )
  11.         for ( k = 0; k < A->n; k++ )
  12.         {
  13.         get_col(A,k,tmp1);
  14.         hhvec(tmp1,k,&beta,tmp1,&(A->me[k][k]));
  15.         hhtrcols(A,k,k+1,tmp1,beta);
  16.         if ( U )
  17.             hhtrcols(U,k,0,tmp1,beta);
  18.         if ( k+1 >= A->n )
  19.             continue;
  20.         get_row(A,k,tmp2);
  21.         hhvec(tmp2,k+1,&beta,tmp2,&(A->me[k][k+1]));
  22.         hhtrrows(A,k+1,k+1,tmp2,beta);
  23.         if ( V )
  24.             hhtrcols(V,k+1,0,tmp2,beta);
  25.         }
  26.     else
  27.         for ( k = 0; k < A->m; k++ )
  28.         {
  29.         get_row(A,k,tmp2);
  30.         hhvec(tmp2,k,&beta,tmp2,&(A->me[k][k]));
  31.         hhtrrows(A,k+1,k,tmp2,beta);
  32.         if ( V )
  33.             hhtrcols(V,k,0,tmp2,beta);
  34.         if ( k+1 >= A->m )
  35.             continue;
  36.         get_col(A,k,tmp1);
  37.         hhvec(tmp1,k+1,&beta,tmp1,&(A->me[k+1][k]));
  38.         hhtrcols(A,k+1,k+1,tmp1,beta);
  39.         if ( U )
  40.             hhtrcols(U,k+1,0,tmp1,beta);
  41.         }
  42.  
  43.     return A;
  44. }
  45.  
  46. /* svd -- returns vector of singular values in d
  47.     -- also updates U and/or V, if one or the other is non-NULL
  48.     -- destroys A */
  49. VEC    *svd(A,U,V,d)
  50. MAT    *A, *U, *V;
  51. VEC    *d;
  52. {
  53.     static VEC    *f=VNULL;
  54.     int    i, limit;
  55.     MAT    *A_tmp;
  56.  
  57.     if ( ! A )
  58.         error(E_NULL,"svd");
  59.     if ( ( U && ( U->m != U->n ) ) || ( V && ( V->m != V->n ) ) )
  60.         error(E_SQUARE,"svd");
  61.     if ( ( U && U->m != A->m ) || ( V && V->m != A->n ) )
  62.         error(E_SIZES,"svd");
  63.  
  64.     A_tmp = m_copy(A,MNULL);
  65.     if ( U != MNULL )
  66.         m_ident(U);
  67.     if ( V != MNULL )
  68.         m_ident(V);
  69.     limit = min(A_tmp->m,A_tmp->n);
  70.     d = v_resize(d,limit);
  71.     f = v_resize(f,limit-1);
  72.     MEM_STAT_REG(f,TYPE_VEC);
  73.  
  74.     bifactor(A_tmp,U,V);
  75.     if ( A_tmp->m >= A_tmp->n )
  76.         for ( i = 0; i < limit; i++ )
  77.         {
  78.         d->ve[i] = A_tmp->me[i][i];
  79.         if ( i+1 < limit )
  80.             f->ve[i] = A_tmp->me[i][i+1];
  81.         }
  82.     else
  83.         for ( i = 0; i < limit; i++ )
  84.         {
  85.         d->ve[i] = A_tmp->me[i][i];
  86.         if ( i+1 < limit )
  87.             f->ve[i] = A_tmp->me[i+1][i];
  88.         }
  89.  
  90.  
  91.     if ( A_tmp->m >= A_tmp->n )
  92.         bisvd(d,f,U,V);
  93.     else
  94.         bisvd(d,f,V,U);
  95.  
  96.     M_FREE(A_tmp);
  97.  
  98.     return d;
  99. }
  100.  
  101. FileDataŵsymmeigOEÿÿÿüz¶(4¿
  102. /**************************************************************************
  103. **
  104. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  105. **
  106. **                 Meschach Library
  107. ** 
  108. ** This Meschach Library is provided "as is" without any express 
  109. ** or implied warranty of any kind with respect to this software. 
  110. ** In particular the authors shall not be liable for any direct, 
  111. ** indirect, special, incidental or consequential damages arising 
  112. ** in any way from use of the software.
  113. ** 
  114. ** Everyone is granted permission to copy, modify and redistribute this
  115. ** Meschach Library, provided:
  116. **  1.  All copies contain this copyright notice.
  117. **  2.  All modified copies shall carry a notice stating who
  118. **      made the last modification and the date of such modification.
  119. **  3.  No charge is made for this software or works derived from it.  
  120. **      This clause shall not be construed as constraining other software
  121. **      distributed on the same medium as this software, nor is a
  122. **      distribution fee considered a charge.
  123. **
  124. ***************************************************************************/
  125.  
  126.  
  127. /*
  128.     File containing routines for symmetric eigenvalue problems
  129. */
  130.  
  131. #include    <stdio.h>
  132. #include    <math.h>
  133. #include    "matrix.h"
  134. #include        "matrix2.h"
  135.  
  136.  
  137. static char rcsid[] = "$Id: symmeig.c,v 1.5 1994/02/16 03:23:39 des Exp $";
  138.  
  139.  
  140.  
  141. #define    SQRT2    1.4142135623730949
  142. #define    sgn(x)    ( (x) >= 0 ? 1 : -1 )
  143.  
  144. /* trieig -- finds eigenvalues of symmetric tridiagonal matrices
  145.     -- matrix represented by a pair of vectors a (diag entries)
  146.         and b (sub- & super-diag entries)
  147.     -- eigenvalues in a on return */
  148. VEC    *trieig(a,b,Q)
  149. VEC    *a, *b;
  150. MAT    *Q;
  151. {
  152.     int    i, i_min, i_max, n, split;
  153.     Real    *a_ve, *b_ve;
  154.     Real    b_sqr, bk, ak1, bk1, ak2, bk2, z;
  155.     Real    c, c2, cs, s, s2, d, mu;
  156.  
  157.     if ( ! a || ! b )
  158.         error(E_NULL,"trieig");
  159.     if ( a->dim != b->dim + 1 || ( Q && Q->m != a->dim ) )
  160.         error(E_SIZES,"trieig");
  161.     if ( Q && Q->m != Q->n )
  162.         error(E_SQUARE,"trieig");
  163.  
  164.     n = a->dim;
  165.     a_ve = a->ve;        b_ve = b->ve;
  166.  
  167.     i_min = 0;
  168.     while ( i_min < n )        /* outer while loop */
  169.     {
  170.         /* find i_max to suit;
  171.             submatrix i_min..i_max should be irreducible */
  172.         i_max = n-1;
  173.         for ( i = i_min; i < n-1; i++ )
  174.             if ( b_ve[i] == 0.0 )
  175.             {    i_max = i;    break;    }
  176.         if ( i_max <= i_min )
  177.         {
  178.             /* printf("# i_min = %d, i_max = %d\n",i_min,i_max); */
  179.             i_min = i_max + 1;
  180.             continue;    /* outer while loop */
  181.         }
  182.  
  183.         /* printf("# i_min = %d, i_max = %d\n",i_min,i_max); */
  184.  
  185.         /* repeatedly perform QR method until matrix splits */
  186.         split = FALSE;
  187.         while ( ! split )        /* inner while loop */
  188.         {
  189.  
  190.             /* find Wilkinson shift */
  191.             d = (a_ve[i_max-1] - a_ve[i_max])/2;
  192.             b_sqr = b_ve[i_max-1]*b_ve[i_max-1];
  193.             mu = a_ve[i_max] - b_sqr/(d + sgn(d)*sqrt(d*d+b_sqr));
  194.             /* printf("# Wilkinson shift = %g\n",mu); */
  195.  
  196.             /* initial Givens' rotation */
  197.             givens(a_ve[i_min]-mu,b_ve[i_min],&c,&s);
  198.             s = -s;
  199.             /* printf("# c = %g, s = %g\n",c,s); */
  200.             if ( fabs(c) < SQRT2 )
  201.             {    c2 = c*c;    s2 = 1-c2;    }
  202.             else
  203.             {    s2 = s*s;    c2 = 1-s2;    }
  204.             cs = c*s;
  205.             ak1 = c2*a_ve[i_min]+s2*a_ve[i_min+1]-2*cs*b_ve[i_min];
  206.             bk1 = cs*(a_ve[i_min]-a_ve[i_min+1]) +
  207.                         (c2-s2)*b_ve[i_min];
  208.             ak2 = s2*a_ve[i_min]+c2*a_ve[i_min+1]+2*cs*b_ve[i_min];
  209.             bk2 = ( i_min < i_max-1 ) ? c*b_ve[i_min+1] : 0.0;
  210.             z  = ( i_min < i_max-1 ) ? -s*b_ve[i_min+1] : 0.0;
  211.             a_ve[i_min] = ak1;
  212.             a_ve[i_min+1] = ak2;
  213.             b_ve[i_min] = bk1;
  214.             if ( i_min < i_max-1 )
  215.             b_ve[i_min+1] = bk2;
  216.             if ( Q )
  217.             rot_cols(Q,i_min,i_min+1,c,-s,Q);
  218.             /* printf("# z = %g\n",z); */
  219.             /* printf("# a [temp1] =\n");    v_output(a); */
  220.             /* printf("# b [temp1] =\n");    v_output(b); */
  221.  
  222.             for ( i = i_min+1; i < i_max; i++ )
  223.             {
  224.             /* get Givens' rotation for sub-block -- k == i-1 */
  225.             givens(b_ve[i-1],z,&c,&s);
  226.             s = -s;
  227.             /* printf("# c = %g, s = %g\n",c,s); */
  228.  
  229.             /* perform Givens' rotation on sub-block */
  230.                 if ( fabs(c) < SQRT2 )
  231.                 {    c2 = c*c;    s2 = 1-c2;    }
  232.                 else
  233.                 {    s2 = s*s;    c2 = 1-s2;    }
  234.                 cs = c*s;
  235.             bk  = c*b_ve[i-1] - s*z;
  236.             ak1 = c2*a_ve[i]+s2*a_ve[i+1]-2*cs*b_ve[i];
  237.             bk1 = cs*(a_ve[i]-a_ve[i+1]) +
  238.                         (c2-s2)*b_ve[i];
  239.             ak2 = s2*a_ve[i]+c2*a_ve[i+1]+2*cs*b_ve[i];
  240.             bk2 = ( i+1 < i_max ) ? c*b_ve[i+1] : 0.0;
  241.             z  = ( i+1 < i_max ) ? -s*b_ve[i+1] : 0.0;
  242.             a_ve[i] = ak1;    a_ve[i+1] = ak2;
  243.             b_ve[i] = bk1;
  244.             if ( i < i_max-1 )
  245.                 b_ve[i+1] = bk2;
  246.             if ( i > i_min )
  247.                 b_ve[i-1] = bk;
  248.             if ( Q )
  249.                 rot_cols(Q,i,i+1,c,-s,Q);
  250.                 /* printf("# a [temp2] =\n");    v_output(a); */
  251.                 /* printf("# b [temp2] =\n");    v_output(b); */
  252.             }
  253.  
  254.             /* test to see if matrix should be split */
  255.             for ( i = i_min; i < i_max; i++ )
  256.             if ( fabs(b_ve[i]) < MACHEPS*
  257.                     (fabs(a_ve[i])+fabs(a_ve[i+1])) )
  258.             {   b_ve[i] = 0.0;    split = TRUE;    }
  259.  
  260.             /* printf("# a =\n");    v_output(a); */
  261.             /* printf("# b =\n");    v_output(b); */
  262.         }
  263.     }
  264.  
  265.     return a;
  266. }
  267.  
  268. /* symmeig -- computes eigenvalues of a dense symmetric matrix
  269.     -- A **must** be symmetric on entry
  270.     -- eigenvalues stored in out
  271.     -- Q contains orthogonal matrix of eigenvectors
  272.     -- returns vector of eigenvalues */
  273. VEC    *symmeig(A,Q,out)
  274. MAT    *A, *Q;
  275. VEC    *out;
  276. {
  277.     int    i;
  278.     static MAT    *tmp = MNULL;
  279.     static VEC    *b   = VNULL, *diag = VNULL, *beta = VNULL;
  280.  
  281.     if ( ! A )
  282.         error(E_NULL,"symmeig");
  283.     if ( A->m != A->n )
  284.         error(E_SQUARE,"symmeig");
  285.     if ( ! out || out->dim != A->m )
  286.         out = v_resize(out,A->m);
  287.  
  288.     tmp  = m_copy(A,tmp);
  289.     b    = v_resize(b,A->m - 1);
  290.     diag = v_resize(diag,(u_int)A->m);
  291.