home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 1 / ARM_CLUB_CD.iso / contents / apps / clib / progs / meschach / !Meschach / c / tutadv < prev    next >
Encoding:
Text File  |  1994-05-19  |  4.4 KB  |  180 lines
  1.     if ( i - l > r - i )
  2.         {    stack[++sp] = l;    stack[++sp] = i-1;    l = i+1;    }
  3.         else
  4.         {    stack[++sp] = i+1;  stack[++sp] = r;    r = i-1;    }
  5.     }
  6.     if ( sp < 0 )
  7.         break;
  8.     r = stack[sp--];    l = stack[sp--];
  9.     }
  10. }
  11.  
  12.  
  13. /* bisvd -- svd of a bidiagonal m x n matrix represented by d (diagonal) and
  14.             f (super-diagonals)
  15.     -- returns with d set to the singular values, f zeroed
  16.     -- if U, V non-NULL, the orthogonal operations are accumulated
  17.         in U, V; if U, V == I on entry, then SVD == U^T.A.V
  18.         where A is initial matrix
  19.     -- returns d on exit */
  20. VEC    *bisvd(d,f,U,V)
  21. VEC    *d, *f;
  22. MAT    *U, *V;
  23. {
  24.     int    i, j, n;
  25.     int    i_min, i_max, split;
  26.     Real    c, s, shift, size, z;
  27.     Real    d_tmp, diff, t11, t12, t22, *d_ve, *f_ve;
  28.  
  29.     if ( ! d || ! f )
  30.         error(E_NULL,"bisvd");
  31.     if ( d->dim != f->dim + 1 )
  32.         error(E_SIZES,"bisvd");
  33.     n = d->dim;
  34.     if ( ( U && U->n < n ) || ( V && V->m < n ) )
  35.         error(E_SIZES,"bisvd");
  36.     if ( ( U && U->m != U->n ) || ( V && V->m != V->n ) )
  37.         error(E_SQUARE,"bisvd");
  38.  
  39.  
  40.     if ( n == 1 )
  41.         return d;
  42.     d_ve = d->ve;    f_ve = f->ve;
  43.  
  44.     size = v_norm_inf(d) + v_norm_inf(f);
  45.  
  46.     i_min = 0;
  47.     while ( i_min < n )    /* outer while loop */
  48.     {
  49.         /* find i_max to suit;
  50.         submatrix i_min..i_max should be irreducible */
  51.         i_max = n - 1;
  52.         for ( i = i_min; i < n - 1; i++ )
  53.         if ( d_ve[i] == 0.0 || f_ve[i] == 0.0 )
  54.         {   i_max = i;
  55.             if ( f_ve[i] != 0.0 )
  56.             {
  57.             /* have to ``chase'' f[i] element out of matrix */
  58.             z = f_ve[i];    f_ve[i] = 0.0;
  59.             for ( j = i; j < n-1 && z != 0.0; j++ )
  60.             {
  61.                 givens(d_ve[j+1],z, &c, &s);
  62.                 s = -s;
  63.                 d_ve[j+1] =  c*d_ve[j+1] - s*z;
  64.                 if ( j+1 < n-1 )
  65.                 {
  66.                 z         = s*f_ve[j+1];
  67.                 f_ve[j+1] = c*f_ve[j+1];
  68.                 }
  69.                 if ( U )
  70.                 rot_rows(U,i,j+1,c,s,U);
  71.             }
  72.             }
  73.             break;
  74.         }
  75.         if ( i_max <= i_min )
  76.         {
  77.         i_min = i_max + 1;
  78.         continue;
  79.         }
  80.         /* printf("bisvd: i_min = %d, i_max = %d\n",i_min,i_max); */
  81.  
  82.         split = FALSE;
  83.         while ( ! split )
  84.         {
  85.         /* compute shift */
  86.         t11 = d_ve[i_max-1]*d_ve[i_max-1] +
  87.             (i_max > i_min+1 ? f_ve[i_max-2]*f_ve[i_max-2] : 0.0);
  88.         t12 = d_ve[i_max-1]*f_ve[i_max-1];
  89.         t22 = d_ve[i_max]*d_ve[i_max] + f_ve[i_max-1]*f_ve[i_max-1];
  90.         /* use e-val of [[t11,t12],[t12,t22]] matrix
  91.                 closest to t22 */
  92.         diff = (t11-t22)/2;
  93.         shift = t22 - t12*t12/(diff +
  94.             sgn(diff)*sqrt(diff*diff+t12*t12));
  95.  
  96.         /* initial Givens' rotation */
  97.         givens(d_ve[i_min]*d_ve[i_min]-shift,
  98.             d_ve[i_min]*f_ve[i_min], &c, &s);
  99.  
  100.         /* do initial Givens' rotations */
  101.         d_tmp         = c*d_ve[i_min] + s*f_ve[i_min];
  102.         f_ve[i_min]   = c*f_ve[i_min] - s*d_ve[i_min];
  103.         d_ve[i_min]   = d_tmp;
  104.         z             = s*d_ve[i_min+1];
  105.         d_ve[i_min+1] = c*d_ve[i_min+1];
  106.         if ( V )
  107.             rot_rows(V,i_min,i_min+1,c,s,V);
  108.         /* 2nd Givens' rotation */
  109.         givens(d_ve[i_min],z, &c, &s);
  110.         d_ve[i_min]   = c*d_ve[i_min] + s*z;
  111.         d_tmp         = c*d_ve[i_min+1] - s*f_ve[i_min];
  112.         f_ve[i_min]   = s*d_ve[i_min+1] + c*f_ve[i_min];
  113.         d_ve[i_min+1] = d_tmp;
  114.         if ( i_min+1 < i_max )
  115.         {
  116.             z             = s*f_ve[i_min+1];
  117.             f_ve[i_min+1] = c*f_ve[i_min+1];
  118.         }
  119.         if ( U )
  120.             rot_rows(U,i_min,i_min+1,c,s,U);
  121.  
  122.         for ( i = i_min+1; i < i_max; i++ )
  123.         {
  124.             /* get Givens' rotation for zeroing z */
  125.             givens(f_ve[i-1],z, &c, &s);
  126.             f_ve[i-1] = c*f_ve[i-1] + s*z;
  127.             d_tmp     = c*d_ve[i] + s*f_ve[i];
  128.             f_ve[i]   = c*f_ve[i] - s*d_ve[i];
  129.             d_ve[i]   = d_tmp;
  130.             z         = s*d_ve[i+1];
  131.             d_ve[i+1] = c*d_ve[i+1];
  132.             if ( V )
  133.             rot_rows(V,i,i+1,c,s,V);
  134.             /* get 2nd Givens' rotation */
  135.             givens(d_ve[i],z, &c, &s);
  136.             d_ve[i]   = c*d_ve[i] + s*z;
  137.             d_tmp     = c*d_ve[i+1] - s*f_ve[i];
  138.             f_ve[i]   = c*f_ve[i] + s*d_ve[i+1];
  139.             d_ve[i+1] = d_tmp;
  140.             if ( i+1 < i_max )
  141.             {
  142.             z         = s*f_ve[i+1];
  143.             f_ve[i+1] = c*f_ve[i+1];
  144.             }
  145.             if ( U )
  146.             rot_rows(U,i,i+1,c,s,U);
  147.         }
  148.         /* should matrix be split? */
  149.         for ( i = i_min; i < i_max; i++ )
  150.             if ( fabs(f_ve[i]) <
  151.                 MACHEPS*(fabs(d_ve[i])+fabs(d_ve[i+1])) )
  152.             {
  153.             split = TRUE;
  154.             f_ve[i] = 0.0;
  155.             }
  156.             else if ( fabs(d_ve[i]) < MACHEPS*size )
  157.             {
  158.             split = TRUE;
  159.             d_ve[i] = 0.0;
  160.             }
  161.             /* printf("bisvd: d =\n");    v_output(d); */
  162.             /* printf("bisvd: f = \n");    v_output(f); */
  163.         }
  164.     }
  165.     fixsvd(d,U,V);
  166.  
  167.     return d;
  168. }
  169.  
  170. /* bifactor -- perform preliminary factorisation for bisvd
  171.     -- updates U and/or V, which ever is not NULL */
  172. MAT    *bifactor(A,U,V)
  173. MAT    *A, *U, *V;
  174. {
  175.     int    k;
  176.     static VEC    *tmp1=VNULL, *tmp2=VNULL;
  177.     Real    beta;
  178.  
  179.     if ( ! A )
  180.         e