home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 1 / ARM_CLUB_CD.iso / contents / apps / clib / progs / meschach / !Meschach / c / svd < prev    next >
Encoding:
Text File  |  1994-01-13  |  9.7 KB  |  405 lines
  1. n,i_max); */
  2.             i_min = i_max + 1;
  3.             continue;    /* outer while loop */
  4.         }
  5.  
  6.         /* printf("# i_min = %d, i_max = %d\n",i_min,i_max); */
  7.  
  8.         /* repeatedly perform QR method until matrix splits */
  9.         split = FALSE;
  10.         while ( ! split )        /* inner while loop */
  11.         {
  12.  
  13.             /* find Wilkinson shift */
  14.             d = (a_ve[i_max-1] - a_ve[i_max])/2;
  15.             b_sqr = b_ve[i_max-1]*b_ve[i_max-1];
  16.             mu = a_ve[i_max] - b_sqr/(d + sgn(d)*sqrt(d*d+b_sqr));
  17.             /* printf("# Wilkinson shift = %g\n",mu); */
  18.  
  19.             /* initial Givens' rotation */
  20.             givens(a_ve[i_min]-mu,b_ve[i_min],&c,&s);
  21.             s = -s;
  22.             /* printf("# c = %g, s = %g\n",c,s); */
  23.             if ( fabs(c) < SQRT2 )
  24.             {    c2 = c*c;    s2 = 1-c2;    }
  25.             else
  26.             {    s2 = s*s;    c2 = 1-s2;    }
  27.             cs = c*s;
  28.             ak1 = c2*a_ve[i_min]+s2*a_ve[i_min+1]-2*cs*b_ve[i_min];
  29.             bk1 = cs*(a_ve[i_min]-a_ve[i_min+1]) +
  30.                         (c2-s2)*b_ve[i_min];
  31.             ak2 = s2*a_ve[i_min]+c2*a_ve[i_min+1]+2*cs*b_ve[i_min];
  32.             bk2 = ( i_min < i_max-1 ) ? c*b_ve[i_min+1] : 0.0;
  33.             z  = ( i_min < i_max-1 ) ? -s*b_ve[i_min+1] : 0.0;
  34.             a_ve[i_min] = ak1;
  35.             a_ve[i_min+1] = ak2;
  36.             b_ve[i_min] = bk1;
  37.             if ( i_min < i_max-1 )
  38.             b_ve[i_min+1] = bk2;
  39.             if ( Q )
  40.             rot_cols(Q,i_min,i_min+1,c,-s,Q);
  41.             /* printf("# z = %g\n",z); */
  42.             /* printf("# a [temp1] =\n");    v_output(a); */
  43.             /* printf("# b [temp1] =\n");    v_output(b); */
  44.  
  45.             for ( i = i_min+1; i < i_max; i++ )
  46.             {
  47.             /* get Givens' rotation for sub-block -- k == i-1 */
  48.             givens(b_ve[i-1],z,&c,&s);
  49.             s = -s;
  50.             /* printf("# c = %g, s = %g\n",c,s); */
  51.  
  52.             /* perform Givens' rotation on sub-block */
  53.                 if ( fabs(c) < SQRT2 )
  54.                 {    c2 = c*c;    s2 = 1-c2;    }
  55.                 else
  56.                 {    s2 = s*s;    c2 = 1-s2;    }
  57.                 cs = c*s;
  58.             bk  = c*b_ve[i-1] - s*z;
  59.             ak1 = c2*a_ve[i]+s2*a_ve[i+1]-2*cs*b_ve[i];
  60.             bk1 = cs*(a_ve[i]-a_ve[i+1]) +
  61.                         (c2-s2)*b_ve[i];
  62.             ak2 = s2*a_ve[i]+c2*a_ve[i+1]+2*cs*b_ve[i];
  63.             bk2 = ( i+1 < i_max ) ? c*b_ve[i+1] : 0.0;
  64.             z  = ( i+1 < i_max ) ? -s*b_ve[i+1] : 0.0;
  65.             a_ve[i] = ak1;    a_ve[i+1] = ak2;
  66.             b_ve[i] = bk1;
  67.             if ( i < i_max-1 )
  68.                 b_ve[i+1] = bk2;
  69.             if ( i > i_min )
  70.                 b_ve[i-1] = bk;
  71.             if ( Q )
  72.                 rot_cols(Q,i,i+1,c,-s,Q);
  73.                 /* printf("# a [temp2] =\n");    v_output(a); */
  74.                 /* printf("# b [temp2] =\n");    v_output(b); */
  75.             }
  76.  
  77.             /* test to see if matrix should be split */
  78.             for ( i = i_min; i < i_max; i++ )
  79.             if ( fabs(b_ve[i]) < MACHEPS*
  80.                     (fabs(a_ve[i])+fabs(a_ve[i+1])) )
  81.             {   b_ve[i] = 0.0;    split = TRUE;    }
  82.  
  83.             /* printf("# a =\n");    v_output(a); */
  84.             /* printf("# b =\n");    v_output(b); */
  85.         }
  86.     }
  87.  
  88.     return a;
  89. }
  90.  
  91. /* symmeig -- computes eigenvalues of a dense symmetric matrix
  92.     -- A **must** be symmetric on entry
  93.     -- eigenvalues stored in out
  94.     -- Q contains orthogonal matrix of eigenvectors
  95.     -- returns vector of eigenvalues */
  96. VEC    *symmeig(A,Q,out)
  97. MAT    *A, *Q;
  98. VEC    *out;
  99. {
  100.     int    i;
  101.     static MAT    *tmp = MNULL;
  102.     static VEC    *b   = VNULL, *diag = VNULL, *beta = VNULL;
  103.  
  104.     if ( ! A )
  105.         error(E_NULL,"symmeig");
  106.     if ( A->m != A->n )
  107.         error(E_SQUARE,"symmeig");
  108.     if ( ! out || out->dim != A->m )
  109.         out = v_resize(out,A->m);
  110.  
  111.     tmp  = m_copy(A,tmp);
  112.     b    = v_resize(b,A->m - 1);
  113.     diag = v_resize(diag,(u_int)A->m);
  114.     beta = v_resize(beta,(u_int)A->m);
  115.     MEM_STAT_REG(tmp,TYPE_MAT);
  116.     MEM_STAT_REG(b,TYPE_VEC);
  117.     MEM_STAT_REG(diag,TYPE_VEC);
  118.     MEM_STAT_REG(beta,TYPE_VEC);
  119.  
  120.     Hfactor(tmp,diag,beta);
  121.     if ( Q )
  122.         makeHQ(tmp,diag,beta,Q);
  123.  
  124.     for ( i = 0; i < A->m - 1; i++ )
  125.     {
  126.         out->ve[i] = tmp->me[i][i];
  127.         b->ve[i] = tmp->me[i][i+1];
  128.     }
  129.     out->ve[i] = tmp->me[i][i];
  130.     trieig(out,b,Q);
  131.  
  132.     return out;
  133. }
  134.  
  135. FileDataŵtorture%nEÿÿÿó8öªÜ
  136. /**************************************************************************
  137. **
  138. ** Copyright (C) 1993 - expand row by means of realloc()
  139.    -- type must be TYPE_SPMAT if r is a row of a SPMAT structure,
  140.       otherwise it must be TYPE_SPROW
  141.    -- returns r */
  142. SPROW    *sprow_xpd(r,n,type)
  143. SPROW    *r;
  144. int    n,type;
  145. {
  146.    int    newlen;
  147.    
  148.    if ( ! r ) {
  149.      r = NEW(SPROW);
  150.      if (! r ) 
  151.        error(E_MEM,"sprow_xpd");
  152.      else if ( mem_info_is_on()) {
  153.     if (type != TYPE_SPMAT && type != TYPE_SPROW)
  154.       warning(WARN_WRONG_TYPE,"sprow_xpd");
  155.     mem_bytes(type,0,sizeof(SPROW));
  156.     if (type == TYPE_SPROW)
  157.       mem_numvar(type,1);
  158.      }
  159.    }
  160.  
  161.    if ( ! r->elt )
  162.    {
  163.       r->elt = NEW_A((unsigned)n,row_elt);
  164.       if ( ! r->elt )
  165.     error(E_MEM,"sprow_xpd");
  166.       else if (mem_info_is_on()) {
  167.      mem_bytes(type,0,n*sizeof(row_elt));
  168.       }
  169.       r->len = 0;
  170.       r->maxlen = n;
  171.       return r;
  172.    }
  173.    if ( n <= r->len )
  174.      newlen = max(2*r->len + 1,MINROWLEN);
  175.    else
  176.      newlen = n;
  177.    if ( newlen <= r->maxlen )
  178.    {
  179.       MEM_ZERO((char *)(&(r->elt[r->len])),
  180.            (newlen-r->len)*sizeof(row_elt));
  181.       r->len = newlen;
  182.    }
  183.    else
  184.    {
  185.       if (mem_info_is_on()) {
  186.      mem_bytes(type,r->maxlen*sizeof(row_elt),
  187.              newlen*sizeof(row_elt)); 
  188.       }
  189.       r->elt = RENEW(r->elt,newlen,row_elt);
  190.       if ( ! r->elt )
  191.     error(E_MEM,"sprow_xpd");
  192.       r->maxlen = newlen;
  193.       r->len = newlen;
  194.    }
  195.    
  196.    return r;
  197. }
  198.  
  199. /* sprow_resize -- resize a SPROW variable by means of realloc()
  200.    -- n is a new size
  201.    -- returns r */
  202. SPROW    *sprow_resize(r,n,type)
  203. SPROW    *r;
  204. int    n,type;
  205. {
  206.    if (n < 0)
  207.      error(E_NEG,"sprow_resize");
  208.  
  209.    if ( ! r ) 
  210.      return sprow_get(n);
  211.    
  212.    if (n == r->len)
  213.      return r;
  214.  
  215.    if ( ! r->elt )
  216.    {
  217.       r->elt = NEW_A((unsigned)n,row_elt);
  218.       if ( ! r->elt )
  219.     error(E_MEM,"sprow_resize");
  220.       else if (mem_info_is_on()) {
  221.      mem_bytes(type,0,n*sizeof(row_elt));
  222.       }
  223.       r->maxlen = r->len = n;
  224.       return r;
  225.    }
  226.  
  227.    if ( n <= r->maxlen )
  228.      r->len = n;
  229.    else
  230.    {
  231.       if (mem_info_is_on()) {
  232.      mem_bytes(type,r->maxlen*sizeof(row_elt),
  233.            n*sizeof(row_elt)); 
  234.       }
  235.       r->elt = RENEW(r->elt,n,row_elt);
  236.       if ( ! r->elt )
  237.     error(E_MEM,"sprow_resize");
  238.       r->maxlen = r->len = n;
  239.    }
  240.    
  241.    return r;
  242. }
  243.  
  244.  
  245. /* release a row of a matrix */
  246. int sprow_free(r)
  247. SPROW    *r;
  248. {
  249.    if ( ! r )
  250.      return -1;
  251.  
  252.    if (mem_info_is_on()) {
  253.       mem_bytes(TYPE_SPROW,sizeof(SPROW),0);
  254.       mem_numvar(TYPE_SPROW,-1);
  255.    }
  256.    
  257.    if ( r->elt )
  258.    {
  259.       if (mem_info_is_on()) {
  260.      mem_bytes(TYPE_SPROW,r->maxlen*sizeof(row_elt),0);
  261.       }
  262.       free((char *)r->elt);
  263.    }
  264.    free((char *)r);
  265.    return 0;
  266. }
  267.  
  268.  
  269. /* sprow_merge -- merges r1 and r2 into r_out
  270.    -- cannot be done in-situ
  271.    -- type must be SPMAT or SPROW depending on
  272.       whether r_out is a row of a SPMAT structure
  273.       or a SPROW variable
  274.    -- returns r_out */
  275. SPROW    *sprow_merge(r1,r2,r_out,type)
  276. SPROW    *r1, *r2, *r_out;
  277. int type;
  278. {
  279.    int    idx1, idx2, idx_out, len1, len2, len_out;
  280.    row_elt    *elt1, *elt2, *elt_out;
  281.    
  282.    if ( ! r1 || ! r2 )
  283.      error(E_NULL,"sprow_merge");
  284.    if ( ! r_out )
  285.      r_out = sprow_get(MINROWLEN);
  286.    if ( r1 == r_out || r2 == r_out )
  287.      error(E_INSITU,"sprow_merge");
  288.    
  289.    /* Initialise */
  290.    len1 = r1->len;    len2 = r2->len;    len_out = r_out->maxlen;
  291.    idx1 = idx2 = idx_out = 0;
  292.    elt1 = r1->elt;    elt2 = r2->elt;    elt_out = r_out->elt;
  293.    
  294.    while ( idx1 < len1 || idx2 < len2 )
  295.    {
  296.       if ( idx_out >= len_out )
  297.       {   /* r_out is too small */
  298.      r_out->len = idx_out;
  299.      r_out = sprow_xpd(r_out,0,type);
  300.      len_out = r_out->len;
  301.      elt_out = &(r_out->elt[idx_out]);
  302.       }
  303.       if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
  304.       {
  305.      elt_out->col = elt1->col;
  306.      elt_out->val = elt1->val;
  307.      if ( elt1->col == elt2->col && idx2 < len2 )
  308.      {    elt2++;        idx2++;    }
  309.      elt1++;    idx1++;
  310.       }
  311.       else
  312.       {
  313.      elt_out->col = elt2->col;
  314.      elt_out->val = elt2->val;
  315.      elt2++;    idx2++;
  316.       }
  317.       elt_out++;    idx_out++;
  318.    }
  319.    r_out->len = idx_out;
  320.    
  321.    return r_out;
  322. }
  323.  
  324. /* sprow_copy -- copies r1 and r2 into r_out
  325.    -- cannot be done in-situ
  326.    -- type must be SPMAT or SPROW depending on
  327.       whether r_out is a row of a SPMAT structure
  328.       or a SPROW variable
  329.    -- returns r_out */
  330. SPROW    *sprow_copy(r1,r2,r_out,type)
  331. SPROW    *r1, *r2, *r_out;
  332. int type;
  333. {
  334.    int    idx1, idx2, idx_out, len1, len2, len_out;
  335.    row_elt    *elt1, *elt2, *elt_out;
  336.    
  337.    if ( ! r1 || ! r2 )
  338.      error(E_NULL,"sprow_copy");
  339.    if ( ! r_out )
  340.      r_out = sprow_get(MINROWLEN);
  341.    if ( r1 == r_out || r2 == r_out )
  342.      error(E_INSITU,"sprow_copy");
  343.    
  344.    /* Initialise */
  345.    len1 = r1->len;    len2 = r2->len;    len_out = r_out->maxlen;
  346.    idx1 = idx2 = idx_out = 0;
  347.    elt1 = r1->elt;    elt2 = r2->elt;    elt_out = r_out->elt;
  348.    
  349.    while ( idx1 < len1 || idx2 < len2 )
  350.    {
  351.       while ( idx_out >= len_out )
  352.       {   /* r_out is too small */
  353.      r_out->len = idx_out;
  354.      r_out = sprow_xpd(r_out,0,type);
  355.      len_out = r_out->maxlen;
  356.      elt_out = &(r_out->elt[idx_out]);
  357.       }
  358.       if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
  359.       {
  360.      elt_out->col = elt1->col;
  361.      elt_out->val = elt1->val;
  362.      if ( elt1->col == elt2->col && idx2 < len2 )
  363.      {    elt2++;        idx2++;    }
  364.      elt1++;    idx1++;
  365.       }
  366.       else
  367.       {
  368.      elt_out->col = elt2->col;
  369.      elt_out->val = 0.0;
  370.      elt2++;    idx2++;
  371.       }
  372.       elt_out++;    idx_out++;
  373.    }
  374.    r_out->len = idx_out;
  375.    
  376.    return r_out;
  377. }
  378.  
  379. /* sprow_mltadd -- sets r_out <- r1 + alpha.r2
  380.    -- cannot be in situ
  381.    -- only for columns j0, j0+1, ...
  382.    -- type must be SPMAT or SPROW depending on
  383.       whether r_out is a row of a SPMAT structure
  384.       or a SPROW variable
  385.    -- returns r_out */
  386. SPROW    *sprow_mltadd(r1,r2,alpha,j0,r_out,type)
  387. SPROW    *r1, *r2, *r_out;
  388. double    alpha;
  389. int    j0, type;
  390. {
  391.    int    idx1, idx2, idx_out, len1, len2, len_out;
  392.    row_elt    *elt1, *elt2, *elt_out;
  393.    
  394.    if ( ! r1 || ! r2 )
  395.      error(E_NULL,"sprow_mltadd");
  396.    if ( r1 == r_out || r2 == r_out )
  397.      error(E_INSITU,"sprow_mltadd");
  398.    if ( j0 < 0 )
  399.      error(E_BOUNDS,"sprow_mltadd");
  400.    if ( ! r_out )
  401.      r_out = sprow_get(MINROWLEN);
  402.    
  403.    /* Initialise */
  404.    len1 = r1->len;    len2 = r2->len;    len_out = r_out->maxlen;
  405.    /* idx1 = idx2 = idx_out = 0; *