ARM Club 1/ARM_CLUB_CD.iso/contents/apps/clib/progs/meschach/!Meschach/c/submat

Text File | 1994-01-13 | 4.4 KB | 182 lines

;

    if ( ! A )
        error(E_NULL,"bifactor");
    if ( ( U && ( U->m != U->n ) ) || ( V && ( V->m != V->n ) ) )
        error(E_SQUARE,"bifactor");
    if ( ( U && U->m != A->m ) || ( V && V->m != A->n ) )
        error(E_SIZES,"bifactor");
    tmp1 = v_resize(tmp1,A->m);
    tmp2 = v_resize(tmp2,A->n);
    MEM_STAT_REG(tmp1,TYPE_VEC);
    MEM_STAT_REG(tmp2,TYPE_VEC);

    if ( A->m >= A->n )
        for ( k = 0; k < A->n; k++ )
        {
        get_col(A,k,tmp1);
        hhvec(tmp1,k,&beta,tmp1,&(A->me[k][k]));
        hhtrcols(A,k,k+1,tmp1,beta);
        if ( U )
            hhtrcols(U,k,0,tmp1,beta);
        if ( k+1 >= A->n )
            continue;
        get_row(A,k,tmp2);
        hhvec(tmp2,k+1,&beta,tmp2,&(A->me[k][k+1]));
        hhtrrows(A,k+1,k+1,tmp2,beta);
        if ( V )
            hhtrcols(V,k+1,0,tmp2,beta);
        }
    else
        for ( k = 0; k < A->m; k++ )
        {
        get_row(A,k,tmp2);
        hhvec(tmp2,k,&beta,tmp2,&(A->me[k][k]));
        hhtrrows(A,k+1,k,tmp2,beta);
        if ( V )
            hhtrcols(V,k,0,tmp2,beta);
        if ( k+1 >= A->m )
            continue;
        get_col(A,k,tmp1);
        hhvec(tmp1,k+1,&beta,tmp1,&(A->me[k+1][k]));
        hhtrcols(A,k+1,k+1,tmp1,beta);
        if ( U )
            hhtrcols(U,k+1,0,tmp1,beta);
        }

    return A;
}

/* svd -- returns vector of singular values in d
    -- also updates U and/or V, if one or the other is non-NULL
    -- destroys A */
VEC    *svd(A,U,V,d)
MAT    *A, *U, *V;
VEC    *d;
{
    static VEC    *f=VNULL;
    int    i, limit;
    MAT    *A_tmp;

    if ( ! A )
        error(E_NULL,"svd");
    if ( ( U && ( U->m != U->n ) ) || ( V && ( V->m != V->n ) ) )
        error(E_SQUARE,"svd");
    if ( ( U && U->m != A->m ) || ( V && V->m != A->n ) )
        error(E_SIZES,"svd");

    A_tmp = m_copy(A,MNULL);
    if ( U != MNULL )
        m_ident(U);
    if ( V != MNULL )
        m_ident(V);
    limit = min(A_tmp->m,A_tmp->n);
    d = v_resize(d,limit);
    f = v_resize(f,limit-1);
    MEM_STAT_REG(f,TYPE_VEC);

    bifactor(A_tmp,U,V);
    if ( A_tmp->m >= A_tmp->n )
        for ( i = 0; i < limit; i++ )
        {
        d->ve[i] = A_tmp->me[i][i];
        if ( i+1 < limit )
            f->ve[i] = A_tmp->me[i][i+1];
        }
    else
        for ( i = 0; i < limit; i++ )
        {
        d->ve[i] = A_tmp->me[i][i];
        if ( i+1 < limit )
            f->ve[i] = A_tmp->me[i+1][i];
        }


    if ( A_tmp->m >= A_tmp->n )
        bisvd(d,f,U,V);
    else
        bisvd(d,f,V,U);

    M_FREE(A_tmp);

    return d;
}

FileDataŵsymmeigOEÿÿÿüz¶(4¿
/**************************************************************************
**
** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
**
**                 Meschach Library
**
** This Meschach Library is provided "as is" without any express
** or implied warranty of any kind with respect to this software.
** In particular the authors shall not be liable for any direct,
** indirect, special, incidental or consequential damages arising
** in any way from use of the software.
**
** Everyone is granted permission to copy, modify and redistribute this
** Meschach Library, provided:
**  1.  All copies contain this copyright notice.
**  2.  All modified copies shall carry a notice stating who
**      made the last modification and the date of such modification.
**  3.  No charge is made for this software or works derived from it.
**      This clause shall not be construed as constraining other software
**      distributed on the same medium as this software, nor is a
**      distribution fee considered a charge.
**
***************************************************************************/


/*
    File containing routines for symmetric eigenvalue problems
*/

#include    <stdio.h>
#include    <math.h>
#include    "matrix.h"
#include        "matrix2.h"


static char rcsid[] = "$Id: symmeig.c,v 1.5 1994/02/16 03:23:39 des Exp $";



#define    SQRT2    1.4142135623730949
#define    sgn(x)    ( (x) >= 0 ? 1 : -1 )

/* trieig -- finds eigenvalues of symmetric tridiagonal matrices
    -- matrix represented by a pair of vectors a (diag entries)
        and b (sub- & super-diag entries)
    -- eigenvalues in a on return */
VEC    *trieig(a,b,Q)
VEC    *a, *b;
MAT    *Q;
{
    int    i, i_min, i_max, n, split;
    Real    *a_ve, *b_ve;
    Real    b_sqr, bk, ak1, bk1, ak2, bk2, z;
    Real    c, c2, cs, s, s2, d, mu;

    if ( ! a || ! b )
        error(E_NULL,"trieig");
    if ( a->dim != b->dim + 1 || ( Q && Q->m != a->dim ) )
        error(E_SIZES,"trieig");
    if ( Q && Q->m != Q->n )
        error(E_SQUARE,"trieig");

    n = a->dim;
    a_ve = a->ve;        b_ve = b->ve;

    i_min = 0;
    while ( i_min < n )        /* outer while loop */
    {
        /* find i_max to suit;
            submatrix i_min..i_max should be irreducible */
        i_max = n-1;
        for ( i = i_min; i < n-1; i++ )
            if ( b_ve[i] == 0.0 )
            {    i_max = i;    break;    }
        if ( i_max <= i_min )
        {
            /* prin